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by 
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Two nanosecond-timescale molecular dynamics simulations of acetylcholinesterase (EC 
3.1.1.7) were analyzed: one unliganded, the other complexed with the snake- venom 
toxin fasciculin 2. These simulation trajectories revealed complex fluctuation of the 
active site gorge. In both simulations, we observe a two-peaked probability distribution 
of the gorge width. Comparing the gorge width with the principal components of motion 
showed that collective motions of many scales contribute to the opening behavior of the 
gorge. Covariance and correlation results, as displayed in the novel "porcupine plots", 
identified the motions of the protein backbone as the gorge opens. 

Fasciculin 2 binds to the gorge entrance of acetylcholinesterase with excellent com- 
plementarity and many polar and hydrophobic interactions; it appears to sterically block 
access of ligands to the gorge. When fasciculin is present, the gorge width distribu- 
tion is altered such that the gorge is more likely to be narrow. Moreover, there are 
large increases in the opening of alternative passages, namely the side door and the 
back door. Finally, the catalytic triad arrangement in the acetylcholinesterase active site 
is disrupted with fasciculin bound. These data support that, in addition to the steric 
obstruction seen in the crystal structure, fasciculin may inhibit acetylcholinesterase by 
combined allosteric and dynamical means. 

On a larger scale, the general infrastructure for solving the time- dependent diffu- 
sion equation using the finite element method has been implemented. Simulations of 
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synaptic transmission were performed using simplified rectilinear models of the neu- 
romuscular junction to demonstrate the effects of synaptic geometries and reactivity 
parameters. Observations from models representing synapses in fast- and slow-twitch 
muscles follow the trends of experimental data. One of our models explains the effects 
of geometry in muscular dystrophy; another demonstrates the capability of our infras- 
tructure to simulate complicated realistic models based on electron microscopy data. 
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Introduction 



I.A The synapse 

A synapse is a junction through which a nerve cell, or a neuron, communicates with 
another cell. There are two types of synapses, categorized by the manner of trans- 
mission: electrical synapses and chemical synapses [4-6]. This dissertation presents 
simulations on some molecules and processes in the cholinergic synapses, the chemical 
synapses which employs acetylcholine (ACh) as the neurotransmitter molecule. The 
neuromuscular junctions, synapses that transmits signal from a neuron to a muscle cell, 
are predominately cholinergic; cholinergic synapses are also present in the central ner- 
vous system. 

I.A.l The chemistry of synaptic transmission 

Signaling from nerve to muscle is mediated by ACh in neuromuscular junctions. The 
mechanism of this process is as follows (Figure 1.1). ACh is synthesized from choline 
and acetyl-CoA and packaged into vesicles in the presynaptic terminal. An action po- 
tential arriving at the terminal of the presynaptic neuron triggers the opening of voltage- 



I.A.l is reproduced in part with permission from Shen, T.; Tai, K.; Henchman, R.; McCammon, J. A. Ace. 
Chem. Res. 2002, 35, 332-340. Copyright © 2002 American Chemical Society. 
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Figure I.l: Schematic of the neuromuscular synapse. (1) Depolarizing of the cell opens 
Ca^+ channels and Ca^"'" enters. (2) Ca^+ stimulates vesicles to release acetylcholine 
(ACh) into the synapse. (3) ACh diffuses and binds to the ACh receptors, opening their 
ion channels. (4) Na+ and K+ diffuse through the channel, depolarizing the postsy- 
naptic terminal. (5) ACh diffuses to acetylcholinesterase (AChE) and is hydrolyzed to 
choline (Ch) and acetate (A). (6) Ch is taken up for resynthesis of ACh and packaged 
into vesicles. 



gated Ca"+ channels. Ca""'" entering the cytosol stimulates the vesicles to release their 
ACh by exocytosis into the synaptic cleft. ACh diffuses across the cleft and binds to 
and activates the postsynaptic receptors, allowing Na+ and K"'" to diffuse across the 
membrane and depolarize the muscle cell. The enzyme acetylcholine (AChE) quickly 
degrades ACh to acetate and choline, deactivating the receptor and tapering off the trans- 
mitted nerve impulse. Choline is taken up again by the presynaptic terminal for resyn- 
thesis of ACh. 

Several diseases are related to dysfunctions of the synapse. For example, one symp- 



torn of Alzheimer's disease is the reduction of ACh concentration; inpatients with myas- 
thenia gravis (serious muscle weakness) there are fewer receptors than normal. In these 
cases, moderate inhibition of AChE is effective in the treatment of the disease to prolong 
the effect of ACh on the receptor. Some AChE inhibitors for treatment of Alzheimer's 
disease are tacrine, rivastigmine, and galantamine [7]. Inhibitors for treatment of myas- 
thenia gravis include pyridostigmine and neostigmine [8]. However, overwhelming in- 
hibition of AChE, particularly by covalent bonding to the active site serine, is invariably 
lethal. Hence, AChE is a prime target for naturally occurring toxins such as the snake 
venom fasciculin 2 (Fas2), pesticides such as parathion and malathion, and chemical 
warfare agents such as sarin, tabun, and VX [9]. 

The work in Chapters II and III explores the dynamics of the AChE molecule and 
mechanism of AChE inhibition by Fas2 using simulation. 

LA.2 The geometries of the neuromuscular synapse 

The most striking feature of the neuromuscular junction, as observed by electron mi- 
croscopy [10, 11] and more recently by electron tomography, is the postsynaptic folds, 
consisting of long valleys and crests directly opposing those presynaptic "active zones" 
(Figure 1.2). In contrast, the presynaptic membrane is smoother in geometry; the ac- 
tive zones therein are arranged in rows, from which the vesicles of neurotransmitters are 
released into the synaptic cleft, the space between the pre- and postsynaptic membranes. 
Not all synapses share the same geometry: synapses of the slow-twitch and the 
fast-twitch muscles in the same organism differ in their geometric structures. Further, 
neuromuscular junctions in dystrophic animals are not geometrically similar to those in 
normal ones. The work in Chapter IV builds a finite element model of the synapse, with 
the geometry as a part of the parameter set, and subsequently seeks to explain the effects 
of differing geometries on synaptic transmission. 




Figure 1.2: Electron micrograph of synaptic junction from mouse gastrocnemius muscle. 



I.B Molecular dynamics 

Molecular dynamics (MD) solves Newton's equation of motion on an atomistic model 
of a molecule to obtain the trajectory of its motion. The positions r(^) of the atoms in 
the molecule are recorded at selected time points t through the simulation; a force field 
is used to describe the interaction between the atoms. 

I.B.l Force field 

The following MD simulations use the AMBER (Assisted Model Building with Energy 
Refinement) force field [12] as the model for biomolecules, and the extended simple 
point charge (SPC/E) model [13] for water molecules. In the AMBER force field, for a 
conformation of the molecule r, the potential energy is modeled as 



£(r) = Y. Kr{r-r,^f+ Y. ^^{9 - 8,^)' 

bonds angles 



v„ 



Y ^[l+cos(«<^-7)] + ^ 



dihedrals 






IJ 
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(1.1) 



which takes into account the contributions from bond lengths and angles, dihedral an- 
gles, van der Waals interactions, and electrostatic interactions. The analytical form of 
the force on the molecule can then be derived using the fact that it is the negative gradient 
of the potential energy, 

-dE{r) 



dv ' 
The acceleration a is obtained using Newton's equation of motion 



(1.2) 



F = ma. (1.3) 

LB.2 Time-stepping 

The position and acceleration at a certain time / are used in the Verlet leap-frog method 
[14-16] to calculate the atom positions after one time-step dt: 

r{t-\-dT) = r(0 + 5?v(/ + i50 (1.4) 

where the velocity at time / + 2^/ is 

v(f+i-50 = \{t-\5t) + 5ta{t) (1.5) 

which requires the velocity at time t — \5t. For calculating the kinetic energy, tempera- 
ture and other properties which requires the velocity at time t, the expression for current 
velocity is used: 

v(0 = ^(v(/ + i50 + v(/-^50). (1.6) 

LB.3 Constrained dynamics 

It has been shown that constraining the bond lengths between hydrogens and the heavier 
atoms, using the SHAKE algorithm [17, 18] described below, can speed up the time- 
stepping calculation while maintaining the dynamical properties observed in the trajec- 
tory on time scales larger than 0.05 ps. 



Used after each unconstrained time-stepping is performed, the SHAKE algorithm 
iterates through a Hst of constraints starting from the first, and adjusts the positions of 
the atoms as each constraint is being applied. Thus, as the ^th constraint, which fixes the 
bond length between atoms / and j, is about to be applied, any previous ^'th constraint 
with k' <k has already been applied; we call this the previous atom position. If there 
is no previous constraint, the previous atom position is the position resulting from the 
unconstrained time-stepping. 

Let nil be the mass of atom I = i^j\ d be the i-j target bond length for this ^th 
constraint; r be the distance between atoms / and j before the last unconstrained time 
step, which satisfies the constraint; r' be the distance after all ^'th constraints with k' <k 
have been applied. The SHAKE algorithm make the following adjustment to the previous 
atom position for atom l: 

5yi = — , l = ij (1.7) 

m 

where 

(1.8) 



2(^ + ^lr-r^ 



It is likely that by making this adjustment, some of the previous constraints will become 
broken. The SHAKE algorithm is repeated until all constraints are satisfied within a 
specified tolerance, usually 10~^ of the target bond length. 

LB.4 Constant temperature and pressure 

To perform M D simulations under the condition of constant temperature and constant 
pressure, the system is coupled to an external bath. The extent of such coupling is de- 
termined by a time constant. The following M D simulations use an external temperature 
and pressure bath that is weakly coupled to the system. 

During each time-step, the velocity of the atoms is proportionally scaled such that 
the change in temperature over this time-step 4^ satisfies 

dT AT 1 , , 

dr-57 ^ ^(^0-^) (1-9) 



where T is the system temperature, Tq is the temperature of the bath, and Ty is the 
temperature-coupling time constant. Recall that the velocity of the atoms v (each with 
mass m) and the system temperature are related through the kinetic energy, namely. 



atoms 



1 2 
—mv 

2 



3 ^ 

2 ^ 



(1.10) 



where A'^ is the total number of atoms, and k^ the Boltzmann constant. It follows that 
the velocity of each atom should be scaled from y to Av, where 



A 



i+^r^-i 



Tj \T 



2 



(1.11) 



Similarly, isotropic pressure coupling is achieved by scaling the boundary sizes and 
the coordinates of the atoms from / to fil, with 



M 



^- — {Po-P) 



(1.12) 



where p is the system pressure, p^ the pressure of the external bath, and Tp the pressure- 
coupling time constant. 



LB.5 Electrostatics with periodic boundary condition 

When periodic boundary conditions are used in a molecular dynamics simulation, van 
der Waals and other short-range interactions can be easily calculated within the unit 
box. However, electrostatic interactions are long-range and usually reach beyond the 
neighboring boxes. To calculate the interaction between an ion and its periodic images, 
the Ewald sum [15] is used. 

For A'^ ions in a cubic unit box of edge length L in a lattice, the electrostatic potential 
energy among all ions in all copies of the unit box is 



A' A' 



r 



II 



1 / i\ i\ 



87re, 



z-z ■ 



r,.. + n 



-1 



(1.13) 



where z- and Zj are the charges on the ion, and the sum over n is the sum over all lattice 
points n = (njL,nyL,n^L), with rix^ny^n^ being integers; the prime excludes summing 



8 



over the cases with / = j for n = 0. For faster convergence in the calculation of this 
potential energy, the Ewald summation is used instead: 



1 Ji ^ / ^ ' erfc(K-|r;; + n|) 

^zz _ ——Y y y z-z- — -—^ — - 

87re.,ee. l,.f-/' ^ r,... + n 



-0/=ly=l V|n|=0 1^7 
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1 «^/_ 1 /-k' 
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II 
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(1.14) 



III 

In Equation 1.14, term I sums, in real space, the original charges with screening charges 
that are by construction Gaussian; term II cancels the contribution from the screening 
charges; term III corrects the self-energy of the screening charges that is included in 
term I. The value K is an arbitrary parameter chosen to optimize the summation. In 
implementation, particle-mesh Ewald [19], an algorithm that scales as NlogN, is used. 

LB.6 Implementation: parallel supercomputing 

The following MD simulations have been carried out using NWChem [20] on two su- 
percomputers, the Cray T3E and the "Blue Horizon". The Blue Horizon is an IBM 
Scalable POWERparallel (SP) machine. Each of the supercomputers consists of multiple 
processors connected by switches for inter-processor communication. 

NWChem calculates the MD trajectory using multiple processors in parallel. The 
calculation work-load is distributed to each of the processors, and the results are then 
collected by a processor and written out to disk storage. The work-load is distributed 
using spatial domain decomposition: the molecular system on which MD is to be per- 
formed is spatially divided into smaller boxes. Each processor is responsible for a box. 
The length of an M D time step is the time taken by the processor that last finishes the 
calculation, plus the time for the inter-processor communication during the distribution 



of work-load and collection of results. The speed or efficiency of the MD simulation can 
be improved by load- balancing schemes. 

I.C Finite element method 

The finite element method [21] solves partial differential equations, of which the diffu- 
sion equation is an instance. 

I.C.l Diffusion equation 

The diffusion equation is of the form 

VDVc = f in QeR"^ (1.15) 

where c, a function of position, is the concentration of the species of interest; D is the 
diffusion constant; and R is the set of all real numbers. Q. is called the domain; d is then 
the number of dimensions of the domain. A boundary condition on the boundary dQ. 
needs to be specified. The Dirichlet boundary condition is 



c 



D 



on do.. (1.16) 



where c^ can be a function of position or time. The Neumann, or refiective, boundary 
condition is 

n-Vc = on dQ.. (1.17) 

where n is the exterior unit normal vector of the boundary. It is admissible to have 
parts of the boundary be Neumann, and the other part be Dirichlet. In the steady-state 
diffusion problem, we simply have f = 0; in the time- dependent problem, we have 
/ = ^, where t is time. 

LC.2 Weak formulation 

Consider the one-dimensional model diffusion problem 

V-Z)Vc(x) = / in ^=[0,1] (1.18) 



10 

with the Dirichlet boundary condition 

c{x) = on ^^ = {0,1}. (1.19) 

It is beneficial to transform this into a weak formulation by introducing a test function 
V : ^ 1-^ R with v = on dQ.. Then, we multiply both sides of (1.18) by v and integrate 
over Q.. 



/ {V-DVc)vdx = fvdx. 

■m Jo. 



The divergence theorem says 



We then define 



A{c,v) := / DVc-Vvdx 
Jo. 

F{v) := [fvdx 
■Ja 



(1.20) 



/ {V-w)zdx = - / w-Vzdx+ l wz-fidS. (1.21) 

Ja Ja Jda 

Setting w = DVc and z = v, we have 

/ {V-DVc)vdx= - / Z)Vc-Vvdx+ J) (DVc)vhdS. (1.22) 

Ja Ja Jda 

Since we have chosen 

y = on do., (1.23) 

we have 

/ {DVc)vhdS = 0, (1.24) 

Jda 

and (1.22) becomes simply 

/ (V • DVc)vdx =- f DVc • Vv dx. (1.25) 

Ja Ja 

Plugging this into (1.20), we have 

- I DVc-Vvdx = I fvdx. (1.26) 

Ja Ja 



(1.27) 
(1.28) 
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to get the weak formulation 



-A{c,v) = F{v). 



(1.29) 



It can be proven that any solution that satisfies the model problem (1.18) also satis- 
fies the following weak-form problem: Find u ^ Hq (Q) such that the weak formulation 
(1.29) is satisfied Vv ^ Hi , where the definitions of the Hilbert spaces are 



H' 



I//1 



H^{0.) := IceH^ c = at x = 0,x=l| 



and the norms and semi-norm are defined as 



^1 ^dc\^ 
— I dx 
dx 



^ 1 



IL^ 



I//1 



^1 



c dx 



I^IIl2 + I^I//i- 



(1.30) 
(1.31) 

(1.32) 

(1.33) 
(1.34) 



LC.3 Galerkin approximation 

The Galerkin method finds an approximation of the true solution in a subset V^^ of the 
Hilbert space (Hq = H^ {Q.) for the model problem). The problem now becomes: Find 
Cj^ e Yf^ C H^ (Q.) such that 

-A(c„vJ=F(v,),Vv^eV^. (1.35) 

V^^ is constructed as a span of a set of bases or "elements"; that is, 

V^ = span{(^^ , (^2 ' • ■ • ' */*«} (1-36) 

These bases (j)- are defined using the following procedure: First, the domain Q. = [0, 1] 
is tesselated into a uniform mesh, with vertices x- = ih where h = 1/N and i = 0,. . . ,N. 
Then, the elements (j)- is thus defined: 

^(x — x^-_l j, X E |.-^,-_] 5-^J 
<^ii^)=l )^(x,-+|-x), xe(x,,x,.^J (1.37) 

otherwise 
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These kinds of elements have three important properties: (a) local support, (b) continu- 
ous, and (c) within the interval between two neighboring vertices, the element is simply 
a linear polynomial. 



Setting 



N 



""h- La ^-V) 



V 



k 



^c.(^^., (L38) 

(/»., /=1,...,7V (1.39) 



in (1.35), we have a set of equations 

-a(|;c^.(^^,(^.J=F((^,.), /=1,...,^ (1.40) 

^-^c/((^.,(^.) = F((^.), /=1,...,7V (1.41) 

This is equivalent to a matrix equation to solve: 

-Ac = F (1.42) 

where each of the entries AJ = A((^',(^-) in the matrix A, each of the entries c- in the 
vector c are to be solved, and each of the entries F^ = F{^-^ in the vector F. Thus, 
a partial differential equation is transformed into a matrix equation which is tractable 
using linear algebra and numerical methods. The solution, which is an approximation 
of the concentration, can be recovered using (1.38). 

LC.4 Generalization and approximation theory 

This treatment by the Galerkin method can be generalized to solve partial differential 
equations in higher dimensions and time-dependent problems. The tesselation can be 
non-uniform, and elements can be designed in a variety of manners. The quality of the 
approximation can be evaluated a posteriori and even a priori. 
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Accounts of Chemical Research [1], of which I am a coauthor. 



II 

Analysis of a 10 ns molecular dynamics 
simulation of mouse 
acetylcholinesterase 



... la dimensione del tempo e andata in frantumi, non possiamo vivere o 
pensare se non spezzoni di tempo che s'allontanano ognuno lungo una sua 
traiettoria e subito spariscono. [22] 



II.A Abstract 

A 10 ns molecular dynamics simulation of mouse acetylcholinesterase was analyzed, 
with special attention paid to the fluctuation in the width of the gorge, and opening 
events of the back door. The trajectory was first verified to ensure its stability. We 
defined the gorge proper radius as the measure for the extent of gorge opening. We 
developed an expression of an inter-atom distance representative of the gorge proper 
radius in terms of projections on the principal components. This revealed the fact that 
collective motions of many scales contribute to the opening behavior of the gorge. Co- 



This chapter is reproduced with permission from Tai, K.; Shen, T.; Borjesson, U.; PhiUppopoulos, M.; 
McCammon, J. A. Biophys. J. 2001, 81, 715-724. Copyright © 2001 Biophysical Society. 



14 



15 



variance and correlation results identified the motions of the protein backbone as the 
gorge opens. In the back door region, side chain dihedral angles that define the opening 
were identified. 

II.B Introduction 

Acetylcholinesterase (AChE) is the enzyme responsible for the termination of signaling 
in cholinergic synapses (such as the neuromuscular junction) by degrading the neuro- 
transmitter acetylchohne. AChE has a gorge, 2 nm deep, leading to the catalytic site. 
Molecular dynamics (MD) simulations [14,23] have shown breathing motions of this 
gorge [24-26]. They also showed that an alternative portal providing access to the cat- 
alytic site is present in AChE. This back door, as it is called, may facilitate rapid solvent 
and product removal [24, 27]. 

In this paper, we collect a 10 ns trajectory of mouse AChE (mAChE) using MD 
simulation. We define and calculate the time series for the proper radius of the gorge 
and for the back door opening events, as observed in our 10 ns M D trajectory. We also 
perform principal component analysis for the trajectory. Then we look for correlations 
among the results from our analyses with a view to finding the important collective 
motions in AChE responsible for its opening behaviors. 

II.C Methods 

II.C.l Molecular dynamics simulation 

The previous 1 ns M D simulation of mouse AChE [27] was extended to afford a trajec- 
tory of 10.8 ns. The set-up procedure has been described in the previous paper, and is 
summarized below. The crystal structure of mouse AChE (mAChE) in complex with 
fasciculin 2 (Fas2) [28] (Protein Data Bank identification code: IMAH) was used to 
model unliganded mAChE. Fas2 was removed from the structure and seven residues in 
mAChE that are missing in IM AH (residues 258-264) were repaired by modeling with 
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Insight II (Accelrys, San Diego, California) to give the initial structure. The protein was 
solvated in a cubic box (9.6 nm edges) of pre-equilibrated water molecules. In order to 
neutralize the — 10 charge of mAChE, 9 sodium ions were placed in the solvent at about 
0.5 nm from the protein surface, and one sodium ion was placed in the active site near 
Trp 86 and His 447. The simulation system had a total of 8 289 solute atoms and 75 615 
solvent atoms. 

The simulation was performed in the isothermal-isobaric ensemble. The solvent 
and solute were separately coupled to temperature reservoirs of 298. 15 K with coupling 
times of 0.1 ps. Pressure was restrained to 1 atmosphere with a coupling time of 0.4 ps. 
All minimization and MD simulation steps were performed using NWChem Version 3.0 
[20] with the AMBER 94 force field [12]. Long-range electrostatic interactions were cal- 
culated using particle-mesh Ewald summation [19]. Bond lengths between hydrogens 
and heavy atoms were constrained using SHAKE [17]. The simulation was performed 
on 32 processors of a Cray t3e parallel supercomputer at the San Diego Supercomputer 
Center over a period of 3 years, consuming a total of about 200 processor-months of 
supercomputer time. Frames were collected at 1 ps intervals for the simulation length 
of 10.8 ns, giving 1.08 x 10"^ frames. The first 700 frames of this trajectory were con- 
sidered the equilibration phase and not used in the main analysis for reasons described 
in the Results section; only the subsequent 10 000 frames (10 ns) were used in the main 
analysis. Thus, the last 300 ps of the previous trajectory [27] is the first 300 ps of the 
present 10 ns trajectory admitted for analysis. 

ILC.2 Proper radius of the gorge 

In order to characterize the degree of the gorge opening with a single variable, we de- 
fined the gorge proper radius p for the conformation of each snapshot as the maximum 
radius of a spherical ligand that can go through the gorge from outside the protein to 
reach the bottom. Equivalently, it is the maximum probe radius that still generates a 
solvent accessible surface with a continuous topology. 

In our algorithm, we first generate the Shrake and Rupley surface with a testing 
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probe sphere of given radius; then we try to determine whether the surface generated by 
the atom Oy in residue Ser 203 (one of the residues in the active site; the bottom of the 
gorge) is topologically continuous with the surface of bottleneck region (represented by 
the atoms Leu 76 C^,, Trp 286 Co, andTyr72 O^^) for that probe radius. Using a binary 
search algorithm to decide what will be the next probe radius, we can determine p with 
desired precision. We started with a test value of 160 pm, and assumed p to be bounded 
by 80 pm and 240 pm. With 6 iterations in the binary search algorithm, we achieved 
a final discretization of 5 pm. We calculated p for each snapshot, and the results were 
collected as a time series p{t). 

ILC.3 Back door opening events 

In addition to the gorge, MD results from this and previous [27] simulations of mAChE 
also showed an alternative opening occasionally large enough for at least water molecules 
to pass. This opening, named the back door [24], is conjectured to assist in releasing 
solvent or reaction products. The back door is formed by the residues Trp 86, Gly 448, 
Tyr 449, and He 451. Instead of using multiple probe radii to calculate the proper radius, 
as in the case of the gorge, we used a single probe radius of 140 pm to test whether it 
is open or closed. We always blocked the gorge entrance while probing for back door 
opening events. (Similarly, we blocked the back door region in the gorge calculation 
described above.) Each frame that the algorithm reported to have a back door opening 
was verified visually using Insight 11. The result is expressed as a time series 

{0 if closed; 
(ii.i) 
1 if open. 

ILC.4 Principal components 

Principal component analysis. Considering only the a-carbon atoms, the A'^-residue 
trajectory can be considered as a vector function of time ?, namely 



r(0= r^^(t),r^y{t),...,r^^(t) 



-\T 



(11.2) 
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of size / = 3N, containing the Cartesian coordinates at time t for residue 1 in the x- 
direction, residue 1 in the j-direction, . . . , up to residue A'^ in the z-direction. 

The ij-entry C- ■ of the covariance matrix C is the covariance of the positions for 
two degrees of freedom / and 7, namely, C-- = ((/■,(?) — ('■,)f) yii^) ~ (''/)?)). where 
{•)[ is the time average over the whole trajectory. Principal component analysis (PCA) 
[25, 29-31] diagonalizes C by solving A = T^CT, so as to obtain the diagonal matrix A 
with the diagonal entries being the eigenvalues ranked by magnitude. The cth column 



Vp'^2'---''^/ 



of the transformation matrix T is the cth eigenvector v^, that is, T 

In our analysis, the algorithm for PCA of the Cartesian coordinates of the a-carbon 
atoms was implemented in the Java programming language, with the JAMA matrix pack- 
age (The MathWorks, Natick, Massachusetts, and National Institute of Standards and 
Technology, Gaithersburg, Maryland). 

The first 5 residues in the M D simulation (residues 4 to 8) and the last 5 (residues 
539 to 543) were removed before the PC A in order to avoid terminal motions excessively 
dominating the analysis. Therefore, each frame contained the Cartesian coordinates for 
A'^ = 530 a-carbon atoms (Leu 9 to Lys 538), or / = 3A'^ = 1590 degrees of freedom. 
Each frame in the 10 ns trajectory is superimposed to the first frame using quaternion 
fitting before being admitted into the PCA. This removes the motions in the rotational 
and translational degrees of freedom. 

Using the transformation matrix T, we can convert between the projections along 
the principal components and the Cartesian coordinates using Ar(?) := r(?) — (r)^ = 



tT 



is a vector of size /; the cth entry thereof 



Tp(0, where p(0 = /;, (0,/?2(0. • • • .-P/(0 

is the projection along the principal component c. 

"Backprojecting": expressing a Cc(-Cc( distance in terms of projections along prin- 
cipal components. Knowing the rules for matrix multiplication, we can solve for each 
of the entries in Ar. For example, the deviation from average for residue 124 in the x- 

direction is Ar^^Ax = ^[,i24xPi + ^2,\24xP2 "^ ^ ^f,\24xPf where v- j24^ is the (scalar) 

entry in the /th eigenvector that corresponds to residue 124 in the x-direction, etc. (For 
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simplicity in notation, we omit explicitly writing out that the rs and ps are indeed func- 
tions of time.) 

The squared distance between residues 124 and 338 (see the Results section for the 
reason of this choice of residues) can then be written as 



^124,338(^) 



V I24x ''338X/ + V 124>' ''338>'j + V 124^ ''338^) (H-^) 

((''l2Jf-(''338^)r+^'"l24x-^''338J^ 
+ (('■l24>')r- ('■338y)f +^''l24^-^''338y) 
+ (('■l24^)r- ('"338^)r+A'-124^-^''338^) ' 

We define d-^^ := v, ^24^ - ^,-^338^ and ^^ := {r^^J, - {r^^Jt^ and similarly for >' and z; 
then we can write 

/ 

'■l24.r-''338x = 'S-^ + ^l,x/^l + ^2,^/^2 +••• + ^/,x/^/ = ^-v + L ^^.-^^^ 



c=l 



and similarly for y and z. The square of these values is 

/ / / 



,) = 'Si + 2 ^ <5x5c,xPc + ^ ^ 5c,,x5c2,x/^Cj/^C2, 



c=\ 



etc. By defining 



Re := 

we can rewrite Equation (11.3) as 

-\2 



^ '•— ^x ~^^y +S^ — \'^124,338/^ 
= 2 {^xSc,x + ^ydc,y + ^z^c,z) 



Oc^,xOc^,x + Oc^jOc^j + dc^^z^c^,Z 



(11.4) 
(11.5) 
(11.6) 



t/, 



24,338* 



/ / 



(0 =s^Y.^cPc+ L Lq 



c=l 



C2=l Cj=l 



t- I L-T x"^ C I Y"^ ^O * 



(11.7) 



We can calculate the time-independent 5, R := [Re], and Q := [Qc^c.] simply by 
looking at the average structure (r)^ and the transformation matrix T. Note that the 

-\2 



square of the distance, 



d 



124,338 



(0 



, is not simply a linear combination of the projec- 



tions Pc{t), but also has cross-terms in the product of two projections pc (t)pc (?). 
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Alternatively, this squared distance can be equated to 



d 



124,338 



(0 



n2 



( ^ 



/ 



a=x,y\z a=x,y,z \c=\ 





"// ^ 


1 


L 


X! ^c,a/^c 




a=x,y,z 


y=\ J 





(11.8) 



Equations (II. 7) and (II. 8) have equivalent parts: the squared average distance S 

2 

d 



124,338 



is Ea'ra' ^he linear part Y^^^cPc is '^Y^a^a (Ec^c a/^c); and the cross-term 



part Ec. Ef, 2c, c,/?c,/^c, isE 



ji,2f cjj 



(Ef4,a/^c)^ 



. Because of this last correspondence, we 



expect to see that the cross-term part will always be larger than zero; this is indeed what 
we see in the Results section. 

II.D Results 



II.D.l Stability of trajectory and mobility of residues 

Figure 11. 1 shows the solute potential energy of the 10.8 ns trajectory. After the first 
700 ps, the solute potential energy ended its decreasing trend and fluctuated around 
75.9 MJ-mol~^ This is one of our reasons for discarding the first 700 ps of trajectory 
and starting the equilibrium statistics thereafter. 

Another justification for discarding the first 700 ps comes from Figure II. 2, where 
we plot the solute potential energy with the gorge proper radius p. The points from the 
first 700 ps cluster are distinct from the rest of the trajectory for their high energies and 
small gorge radii. This suggests that the mAChE molecule relaxes from the conforma- 
tion in the crystal structure IMAH with higher energy and smaller gorge proper radius 
to one that has a wider gorge and is lower in energy. 

The root mean square deviation (rmsd) from the crystal structure through the 10 ns 
trajectory is shown in Figure II.3. The rmsd for the a-carbon atoms kept stable around 
120 pm, and that for all heavy atoms maintained at 170 pm. 
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Figure 11. 1: The solute potential energy for the 10.8 ns trajectory. 
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Figure II. 2: The solute potential energy versus the gorge proper radius p. The red dots 
are the first 700 ps of trajectory. The black dots are the following 10 ns of trajectory, 
with the adjacent frames connected with lines. 
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Figure II. 3: Time dependence of the root mean square deviations (rmsd) for the 10 ns 
mAChE MD trajectory, using the mAChE part in the 1 M A H crystal structure as reference. 
Red, rmsd for alpha-carbon atoms; green, rmsd for heavy atoms. 



The temperature, the total energy, the mass density, and the volume during the 
10 ns trajectory remained stable. The energy components were inspected to ensure the 
stability of the trajectory. 

The isotropic temperature (B) factor can be calculated from the mean square fluc- 
tuation (msf) [14, 23] using the equation 



B 



(msf) . 



(11.9) 



The B factor for each residue from the mAChE part in the IM AH crystal structure and 
that calculated from the msf in the M D simulation are shown in Figure II.4. 

The time-averaged rmsd for each residue, using the mAChE part in the 1 M A H crys- 
tal structure as the reference was also calculated; it showed similar features as the sim- 
ulated B factor (data not shown). Overall we do not observe outstanding fluctuations in 
the gorge region over other parts of the protein. 
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Figure 11.4: The B factor for each residue from the mAChE part in the IMAH crystal 
structure. Black, crystal structure; red, calculated from the msf in the 10 ns MD trajec- 
tory. 

ILD.2 Back door opening events 

In the 10 000 frames analyzed, only 78 had back door openings (Figure II.5). We com- 
pared T(/) with the dihedral angle values in the residues that form the back door, namely, 
Trp 86, Gly 448, Tyr 449, and He 45 1 . The most important ones are shown in Figure II. 6. 



ILD.3 Correlation and covariance analyses 

The 10 ns time series for the gorge proper radius p(?) is shown in Figure II. 7. The 
average of the proper radius is 152 pm. The probability distribution over all observed 
values of p is not Gaussian, as previously described [32]. 

We define the correlation between the ith degree of freedom, r-{t), and the gorge 
proper radius p{t) to be 

< (r,.(0 -('•,■),) (P(0-(P),)), 



^J {{>■:{')- in) f)^{(p{t)-{p)rf)^ 



(11.10) 



We calculated a correlation vector with 3 entries thus defined for the a-carbon atom 
of each residue (Figure 11. 8). Also, we calculated the average correlation vector within 
each of the secondary structure elements (a-helices and j8-strands; Figure IL9). 
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Figure 11. 5: Stereographs of two conformations in the back door region. In each con- 
formation, the outside of the protein is on the top; Trp 86 is on the lower-left, Gly 448 
on the lower-right, and Tyr 449 in the middle. The Connolly surface generated with a 
probe radius of 140 pm is shown. Top: the frame at 291 ps, showing a closed back door. 
In this frame, ;t2 ^^^rp 86 is 114°, I// of Gly 448 is ^6°,X\ of Tyr 449 is -121°, X[ of 
lie 451 is -72°, and Xi of lie 451 is 161°. Bottom: the frame at 8996 ps, showing an 
open back door. In this frame, Xi of Trp 86 is 120°, \\f of Gly 448 is 65°, X{ of Tyr 449 
is-122°, ;t:^ of lie 451 is -64°, and ;t:2 of lie 451 is 163°. 
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Figure II. 6: Selected dihedral angles through the 10 ns that are important to the opening 
of the back door. The red lines at the bottom mark the back door opening events. The 
dots are the dihedral angles: X2 ^^ ^rp 86 (green); ij/ of Gly 448 (blue); Xi of Tyr 449 
(pink); X[ of He 451 (cyan); Xo of He 451 (brown). 



2.5 




0.5 



2000 4000 6000 

time / ps 



8000 



10000 



Figure II. 7: The 10 ns time series for the gorge proper radius, p{t). 
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Figure II. 8: The correlation vector for each residue displayed on the initial structure 
(see text for definition). The vectors point from green to red; the lengths are in arbitrary 
units. The AChE molecule is shown in the conventional orientation: the viewer looks 
down the gorge, with the A'^-terminus on the top of the figure and the C-terminus at the 
lower-left corner. The active site Ser 203 is shown in the space-filling model. 
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Figure II. 9: The correlation vector within each secondary structure element displayed on 
the initial structure (see text for definition). The a-helices are shown as green cylinders; 
the j8-strands are shown in yellow. The vectors point from green to red; the lengths are 
in arbitrary units. The AChE molecule is shown in the conventional orientation: the 
viewer looks down the gorge, with the A'^-terminus on the top of the figure and the C- 
terminus at the lower-left corner. The active site Ser 203 is shown in the space-filling 
model. Helix 14 is marked with numerals '14'. 
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Similarly, the average velocity covariance is defined as 

((r,.(0-r,.(/-lps))(p(0-p(/-lps))), (11.11) 

The average velocity covariance vector is shown in Figure 11.10. 

ILD.4 Principal component analysis and backprojecting 

We found that the distance between Phe 338 C^2 ^^id Tyr 124 O^^ and the gorge proper 
radius are highly correlated, with a correlation coefficient of 0.91 (Figure 11.11); the 
correlation coefficient between the distance Phe 338 C^-Tyr 124 C^ and the gorge 
proper radius is 0.55. 



"l24.338V^) 



n2 



, namely, the linear 



We calculated the time-dependent contributions to 
part Y.c^cPc{t) and the cross-term part Y.c^'Lc2Qc^c,Pc^{t)Pc2{t) (Figure 11.12). Note 
that sometimes the linear part and the cross-term part have comparable magnitudes but 
opposite directions, for example, around 8.1 ns. 

lI.E Discussion 

lI.E.l S factors 

The crystallographic B factors include a variety of contributions (e.g., crystal contacts, 
static disorder in the crystal) in addition to that from the dynamical fluctuations, and 
these additional contributions are expected to be substantial for a structure of modest 
resolution [33]. The crystallographic resolution is 0.32 nm (3.2 A) for the structure 
IMAH [28]. Generally the crystallographic B factors were larger than those observed in 
M D because the simulation sampled far fewer conformations of the protein than crystal- 
lography; the termini, however, exhibited larger B factors in the MD simulation where 
the protein is solvated in water rather than packed in a crystal. 

In the original crystal structure [28], loop II of Fas2 interacts with the peripheral 
anionic site of mAChE, thereby blocking substrate access to the active site; loop I of 
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Figure 11.10: The average velocity covariance vector for each residue displayed on the 
initial structure (see text for definition). The vectors point from green to red; the lengths 
are in arbitrary units. The AChE molecule is shown in the conventional orientation: 
the viewer looks down the gorge, with the A'^-terminus on the top of the figure and the 
C-terminus at the lower-left corner. The active site Ser 203 is shown in the space-filling 
model. 
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Figure 11.11: The gorge proper radius p versus the distance between Phe 338 0^2 ^^^ 
Tyr 124 Op^ . The correlation coefficient is 0.91. 
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Fas2 fits in a crevice near the entrance of the gorge. Our simulation starts from this 
crystal structure, with the inhibitor Fas2 removed. However, we do not see larger B 
factors in MD than in the crystal structure for the residues involved in contact with Fas2; 
this is probably due to the same reasons mentioned in the previous paragraph. Nor were 
any outstanding deviations observed in the Fas2 contact interface in the time-averaged 
rmsd for each residue. 

ILE.2 Access to the active site 

Unlike the previous 0.5 ns M D simulation on the Torpedo californica AChE (TcAChE) 
complexed with tacrine [25], our simulation showed a non-Gaussian probability distri- 
bution of the gorge proper radius as previously described [32]. In addition, our average 
radius, 152 pm, was smaller than that in the tacrine-complexed TcAChE simulation 
(about 190 pm). The time series p{t) almost never reached the 240 pm; for a substrate 
as large as acetylcholine, AChE hardly allowed it any spontaneous access to the active 
site during the 10 ns trajectory. 

The presence of the ligand tacrine in the gorge may be the reason for the larger 
proper radius in the tacrine-complexed TcAChE simulation. Favorable electrostatic in- 
teractions between AChE and its substrate may help in overcoming the difficulty of 
squeezing through a narrow gorge. Note that the catalysis of AChE occurs on a mil- 
lisecond time-scale; as long as the frequency of gorge opening is not so low that the 
substrate diffuses away from the gorge entrance before commitment to catalysis, rarity 
of opening in this nanosecond simulation does not preclude diffusion-controlled kinet- 
ics. 

ILE.3 Back door opening events 

The proposition for the existence of the back door has been based on visual observation 
of the crystal structure [34] and supported by conformations sampled by simulations of 
TcAChE [24, 35]. In our simulation, although only 78 frames out of our 10 ns trajectory 
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had back door opening events, the back door opening observed here aligned sequentially 
with that observed in the MD simulation of TcAChE: Trp 86, Gly 448, and Tyr 449 in 
mAChE correspond respectively to Trp 84, Gly 441, and Tyr 442 in JcAChE. This is 
a remarkable fact: this alternative opening, named the back door, is not limited to a 
single species, but has been observed by molecular dynamics simulations in at least two 
species. 

The selected dihedral angles in the region, as shown, have preferred values when 
the back door opens. For most of the opening events, e.g., around 3 ns, X2 ^^ Trp 86 
acquires higher values around 150° than the usual values around 110°. Back door 
opening is more likely to happen when I// of Gly 448 assumes a value around 45° 
rather than around 80°, and when x^ of Tyr 449 is around —140° rather than —120°. 
Four major types of rotamers for the x^ and x^ of lie 451 are observed: (/) (x^ ^ 
— 60°,;i^2 ~ ±180°), or (—gauche, trans); observed most often in this trajectory, this 
seems to be almost a necessary, but not sufficient, condition for back door opening. (//) 
(Xi ~ 30°,;i^2 ~ 60°), or (^gauche, ^gauche); observed sporadically in the first half of 
the trajectory, it is interspersed between opening events. (///) (Xi ~ 60°,;)^2 ~ ~60°), 
or (^gauche, —gauche); observed around 5 ns and 6 ns, this rotamer does not give any 
opening events, (iv) (x^ ~ — 60°,;i^2 ~ 60°), or (—gauche, -\-gauche); observed after 
7 ns, this rotamer also does not give any opening events. 

It is difficult to be conclusive about preferred dihedral angle values for back door 
opening for the following reasons: We only have a small number of opening events. In 
addition, some of the opening events in the last 2 ns of the simulation have violations of 
preferred dihedral angle values; for example, the two conformations shown in Figure IL5 
have subtly different values for these dihedral angles, but the one at 8996 ps has an open 
back door but that at 291 ps does not. 

Site-directed mutagenesis experiments on human AChE [36] and TcAChE [35] do 
not support significant participation of back door traffic in catalytic activities. On the 
other hand, there are recent experimental evidences from two different methods sup- 
porting the back door hypothesis: product clearance through the back door is implied 
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by the crystal structure of a carbamoylated TcAChE [37]; one inhibitory monoclonal 
antibody of Electrophorus AChE is reported to bind to the region of the back door [38]. 
It is our view that one has to be cautious while attempting to interpret and reconcile 
these data concerning the back door, and take into consideration the time-scale of each 
method of observation. Indeed, even this simulation of 10 ns only covers a fraction of a 
thousandth of a catalytic cycle of AChE. Despite consistently observing conformations 
with an open back door in several MD simulations, we cannot ascertain details of the 
involvement of the back door in AChE catalytic function. 

ILE.4 Correlation and covariance analyses 

We have found that the "porcupine plots" (Figures II. 8, II. 9 and 11.10) are very helpful 
in visualizing the concerted motions between parts of the protein and the functionally 
interesting motion, the gorge opening. The correlation vectors in Figure II. 8 reveal how 
much each residue moves in concert with the gorge. Note that a correlation vector does 
not show how much the residue fluctuates or how flexible it is, but how it moves from 
its average position when the gorge changes in proper radius. The moiety of AChE that 
includes the gorge (the half of the protein that is closer to the viewer in the figures) 
has remarkable concerted motion with the gorge proper radius. The residues in this 
moiety generally have correlation vectors pointing away from the gorge. These residues 
apparently move away from the gorge when the gorge opens. Some residues that are 
in the other moiety (closer to the base of the gorge and further from the viewer) have 
correlation vectors that are generally smaller. The discrimination of these concerted 
motions is even more pronounced in Figure 11.10: the residues that construct the gorge 
itself have the largest average velocity covariance vectors, while those farthest from the 
gorge have the smallest average velocity covariance vectors. We introduce Figure II. 9 as 
an easy way to visualize and interpret how the secondary structure elements contribute 
to the opening of the gorge. For example, helix 14 (Gly 335 to Tyr 341) is close to 
the active site gorge; it has a large vector in this figure pointing away from the gorge, 
which means it moves away when the gorge opens. This helix includes Phe 338, one 
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of the two residues whose separation correlates well with the gorge proper radius, as 
described above. Most of the a-helices in the same moiety as the gorge move away 
from the gorge; those further from the gorge show little concerted motion. The helices 
and two large j3-sheets (B and C) in the other moiety, which is at the bottom of the 
gorge, generally have upward motion towards the gorge as the the latter opens. 

We speculate that inhibitors such as Fas2 may decrease the likelihood of gorge 
opening by restricting groups on the surface of AChE that have large concerted motions, 
in addition to sterically occluding the entrance to the gorge. Other allosteric inhibition 
mechanisms are also possible. These remain to be investigated and confirmed by our 
current work on the MD simulation of the complex of Fas2 and mAChE. We understand 
that in vivo AChE exists in cholinergic synapses as oligomers indirectly anchored on the 
post-synaptic membrane. These attachments to the AChE monomer are not represented 
in this simulation; they may also modulate the opening behavior of the gorge. 

ILE.5 Contributions of the principal components in gorge opening 

Several convincing pieces of evidence support a complex opening behavior of the gorge. 
Starting from a naive guess that one single principal component or just a few, represent- 
ing a large collective motion in the protein, has a dominant contribution in the gorge 
opening behavior, we shall now enumerate evidence to demonstrate the contrary. 

There is a high correlation (0.91) between the Phe 338 Cg2~Tyr 124 O^ distance 
and the gorge proper radius; this distance alone may be used to predict the gorge proper 
radius with high confidence. The corresponding a-carbon distance has a lower correla- 
tion coefficient (0.55). The side chains seem to contribute significantly to the opening of 
the gorge. As these particular side chains project directly toward the gorge and actually 
constitute the bottleneck much of the time, it makes sense that the correlation of the side 
chain distance and the gorge radius will be particularly strong. A predictive model of 
the gorge width with high accuracy will thus have to include the side chains into con- 
sideration; for that purpose, it is inappropriate to only concentrate on large backbone 
motions and ignore those of the side chains. On the other hand, the largest fluctuations 
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in the gorge radius can be expected to reflect movement of backbone segments, and we 
investigate these in the current work. Ligands larger than acetylcholine do bind to the ac- 
tive site, including a number of important inhibitors, and with the backbone correlation 
analysis we are beginning to explore the mechanism of such processes. 

The eigenvalues from the PCA, ranked by magnitude, decrease rapidly: the 15th 
largest eigenvalue is less than one-tenth (0.1) of the largest, and the 113rd largest is 
less than one-hundredth (0.01) of the largest. In spite of this, there is not one principal 
component that has an outstandingly high correlation with p(?); in other words, no one 
principal component dominates the opening of the gorge. The naive guess expects dis- 
missible cross-term parts (Y,c Ec Qc c Pc Pc^) in the backprojecting expression (Equa- 
tion II.7); in addition, it expects the the linear part {Y.c^cPc{t)) to have one dominant 
term with a large coefficient R^ while all other terms have small coefficients. Contrary 
to this expectation, when ranked by magnitude, there was not a dominant minority of 
terms that stand out in the list of the Re coefficients; indeed, none of the RcPc{t) terms 
in the linear part are dominant (data not shown). It is not possible to single out principal 
components that dominate the width of the gorge. 

Furthermore, we saw in Figure 11.12 that a model that ignores the contribution of 
either the cross-term part or the linear part to the gorge opening behavior is not justified. 
Both the linear part and the cross-term part contribute significantly to fluctuations in the 
width of the gorge. In sum, the naive guess that a few principal components dominate 
the gorge is not justified. 

We conjecture that this finding is related to our fractal dynamics model previously 
proposed [32]. In that model, the time-dependent fluctuation of gorge proper radius 
reveals fractal dynamics that lacks a characteristic time-scale over the several orders of 
magnitude being observed. We suggested that this behavior is caused by a hierarchy of 
motions on different time-scales acting together. It is likely that these motions are related 
to our present results, where many motions on different length-scales (as revealed by 
PCA) have a similar effect on the gorge behavior. 

As a side note, we did attempt to perform energy landscape analysis [39-41] by 
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Figure 11.13: The solute potential energy projected onto a histogram of the first two 
principal components. 

projecting of the solute potential energy onto the first two principal components obtained 
from PC A (Figure 11.13). For a protein as large as AChE, our 10 000 frames did not 
provide enough sampling to show a meaningful energy landscape. Comparison between 
the solute potential energy and the gorge proper radius was also inconclusive, likely for 
the same reason. 

II.F Conclusions 



The 10 ns MD simulation of mAChE was equilibrated and stable as shown by various 
properties and energy components, and gave reasonable B factors. 

The small fraction of total number of frames in which we observed an open back 
door severely limits the statistics of our analysis. Back door opening seems to prefer 
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certain dihedral angle values for some of the residues in the region. 

We developed the "porcupine plots" to visualize the concerted motions between 
residues in AChE and the gorge width. The plots show that residues lining the gorge 
have the largest velocity covariance with the gorge, and residues in the same moiety as 
the gorge move concertedly away from the gorge when it opens. This clearly reveals 
the residues involved in gorge opening, and gives the magnitude and direction of such 
involvement. 

We have found that the Phe 338 C^2~Tyr 124 Op^ distance is highly predictive of the 
gorge proper radius; the corresponding a-carbon distance has a moderate correlation. 
We have not investigated side chain motions extensively in this paper, but side chains 
will be important in a predictive model of the gorge opening. 

Principal component analysis showed that the gorge proper radius is determined 
by motions on many different length- scales; this is likely to be related to the fractal 
dynamics model we previously proposed. 
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Mechanism of acetylcholinesterase 
inhibition by fasciculin: 5 ns molecular 
dynamics simulation 




III.A Abstract 

Our previous molecular dynamics simulation (10 ns) of mouse acetylcholinesterase (EC 
3.1.1.7) revealed complex fluctuation of the enzyme active site gorge. Now we report a 
5 ns simulation of acetylcholinesterase complexed with fasciculin 2. Fasciculin 2 binds 
to the gorge entrance of acetylcholinesterase with excellent complementarity and many 
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polar and hydrophobic interactions. In this simulation of the protein-protein complex, 
where fasciculin 2 appears to sterically block access of ligands to the gorge, again we 
observe a two-peaked probability distribution of the gorge width. When fasciculin is 
present, the gorge width distribution is altered such that the gorge is more likely to 
be narrow. Moreover, there are large increases in the opening of alternative passages, 
namely the side door (near Thr 75) and the back door (near Tyr 449). Finally, the 
catalytic triad arrangement in the acetylcholinesterase active site is disrupted with fas- 
ciculin bound. These data support that, in addition to the steric obstruction seen in the 
crystal structure, fasciculin may inhibit acetylcholinesterase by combined allosteric and 
dynamical means. 

III.B Introduction 

Acetylcholinesterase (AChE; EC 3.1.1.7) terminates synaptic transmission at cholinergic 
synapses by catalyzing hydrolysis of the neurotransmitter acetylcholine [4]. A gorge, 
2 nm in depth, leads from the surface of the enzyme to its buried active site. Fas- 
ciculin 2 (Fas2), a peptidic three-finger snake toxin (61 residues) from green mamba 
{Dendroaspis angusticeps) venom, is a picomolar inhibitor of mammalian acetylcholin- 
esterases. 

Much kinetics and structural information for the fasciculin 2-mouse acetylcholin- 
esterase complex (Fas2-mAChE) has been collected in the past decade. Some kinetic 
data suggested that the inhibition may be achieved by steric obstruction: Fas2 capping 
the entrance of the gorge, sterically blocking the entry of the ligand [42]. The crystallo- 
graphic structures of complexes of Fas2 with mouse acetylcholinesterase (mAChE) [28], 
Torpedo californica AChE [43], and human AChE [44] all supported this hypothesis, 
and showed excellent surface complementarity between Fas2 and AChE. The contacts 
between Fas2 and mAChE are in turn verified by mutation studies [45-47]. 

In the crystal structures, allosteric effects seem to be subtle at most. However, this 
has not been fully reconciled with kinetics information, which shows residual catalytic 
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activity for the Fas2-bound enzyme. Indeed, inhibition of AChE by fasciculin, despite 
its high affinity and slow rate of dissociation, is not complete, and the fractional residual 
activity reveals that substrates can still enter the active site in Fas2-bound AChE [45, 
48-51]. With certain AChE mutants and substrates whose catalysis is not diffusion 
limited, Fas2 even acts as an allosteric activator [52], connoting possible influence by 
Fas2 on alternative passages and other interaction in addition to direct steric occlusion. 
Moreover, it has been suggested that Fas2 alters the conformation of AChE in the active 
site [49] and in the Q. loop (Cys 69 to Cys 96) [28, 53-55], particularly that of the Trp 86 
side chain [52]. 

We have previously carried out a 10 ns molecular dynamics (MD) simulation of 
mAChE [2]. This has revealed the complex nature of the gorge fluctuations. Collec- 
tive motions on many scales contribute to the opening behavior of the gorge, and two 
distinct states, one narrow and one wide, were found. Correlation results identified the 
motions of many protein residues as the gorge opens. In particular, the residues within 
the mAChE moiety that includes the gorge apparently move away from the gorge en- 
trance when the gorge opens. The opening of alternative passages to the active site was 
found to be infrequent, since less than one-hundredth of the frames collected showed 
opening of alternative passages. These alternative passages are the back door, bounded 
by residues Trp 86, Tyr 449, and He 451 [24]; and the side door, by residues Thr 75, 
Leu 76, and Thr 83 [25,27]. 

To investigate the structural and dynamical effect of Fas2 binding on mAChE, we 
performed another M D simulation of the mAChE complexed with Fas2 and compare 
mAChE properties with the results from the previous apo-mAChE trajectory. 

III.C Methods 

III.C.l Crystal structure 

The crystal structure of the Fas2-mAChE complex used for this study (Protein Data 
Bank identification code: lKU6) was refined to 0.25 nm (2.5 A) resolution. It con- 
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tains one Fas2 molecule (residues 1-61) bound to one mACliE monomer (residues 3- 
541), one GlcNAc moiety linked to mAChE residue Asn 350, one ethylene glycol and 
189 water molecules. The overall deviations from ideal geometry are 0.010 A for bond 
distances and 1 .76° for bond angles. This structure, though virtually identical to the pre- 
vious complex structure IMAH [28] with an average root mean square deviation (rmsd) 
of 0.57 A for 587 C^ atoms, provides significantly improved accuracy in the positions 
of the main and side chains of Fas2 and mAChE and of those located at the complex in- 
terface, and unambiguously reveals the position of a higher number of water molecules. 
A detailed description of this structure including a report on the crystallization and data 
collection conditions will be published elsewhere. 

IILC.2 Molecular dynamics simulation 

The 0.25 nm crystal structure of Fas2-mAChE described above was used to build the 
initial structure (Figures ELI and III. 2). In addition to the Fas2 and mAChE molecules 
themselves, the 189 crystallographic water molecules present in the structure were left 
intact and considered part of the solvent throughout the following preparation. 

To mend the missing or truncated residues and segments of the Fas2-mAChE crys- 
tal structure lKU6, we used mAChE segment Pro 258-Gly 264 from the mAChE crys- 
tal structure (IMAA) [57]; mAChE residues Lys 496, Ala 542, and Thr 543 from the 
previous Fas2-mAChE crystal structure (IM AH); and Fas2 residue Lys 51 from the un- 
liganded Fas2 crystal structure (IFSC) [58]. The mAChE residue Arg 493, which in the 
structure was modeled as an alanine residue due to side chain mobility, was kept un- 
changed, following our previous simulation of mAChE [2, 27]. In the end, the prepared 
structure has the same mAChE atoms as our previous apo-mAChE M D . These structure 
manipulations were carried out using Insight II (Accelrys, San Diego, California). 

The following preparation, minimization and M D procedures were performed using 
NWChem version 4.0 [20], using the AMBER 94 force field [12] for the solutes (proteins 
and counterions), and SPC/E [13] for the water molecules. 

Hydrogen was added to the structure using the prepare module of NWChem. To 
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Figure III.l: Overview of the Fas2-mAChE complex structure (lKU6). The mAChE 
molecule is labeled every 25 residues; the Fas2 molecule, moved from the mAChE 
for the sake of clarity, has residues 10 (in loop I), 30 (in loop II) and 50 (in loop III) 
labeled. The mAChE atoms that are used to define the blocking sphere (see Methods) 
for each passage are marked with small spheres: gorge (Tyr 72 O^^, Leu 76 C^, , Trp 286 
Co); side door (Thr 75 Co); back door (Tyr 449 6-membered ring). Generated using 
Molscript [56]. 
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Figure III.2: Stereograph of the secondary-structure representation the Fas2-mAChE 
crystal structure (lKU6). The mAChE molecule is gradually colored from blue (A'^- 
terminus; residue 4) to light orange (C-terminus; residue 543). The Fas2 molecule 
is closer to the viewer, colored from dark orange (A'^-terminus; residue 1) to red (C- 
terminus; residue 61). Generated using Molscript [56]. 
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neutralize the —6 charge of the Fas2-mAChE complex, 6 sodium counterions were 
added; in addition to the sodium ion in the active site as in the previous simulation 
of apo-mAChE [2, 27], the position of the other 5 counterions are analogous to those 
farthest from Fas2 in the previous mAChE simulation. The longest dimension of the 
molecule is 8.4 nm at this point. The molecule is then solvated in a cubic box of water 
molecules with length 10.5 nm on each side. 

The MD simulation was performed in the isothermal-isobaric ensemble. The sol- 
vent and solute were separately coupled to temperature reservoirs of 298.15 K with 
coupling times of 0.1 ps; pressure was restrained to 1.025 x 10^ Pa with a coupling time 
of 0.5 ps [59]. The M D simulation time step was 2 fs for both equilibration and produc- 
tion. Long-range electrostatic interactions were calculated using particle-mesh Ewald 
summation [19]. Bond lengths between hydrogens and heavy atoms were constrained 
using SHAKE [17, 18]. 

The water molecules were relaxed using steepest descent for 5 000 steps, and equi- 
librated at 298.15 K for 100 ps. The whole M D simulation system contained 540 residues 
(residues 4 to 543; 8 279 atoms) for mAChE, 61 residues (906 atoms) for Fas2, 6 sodium 
counterions, and 35 796 solvent water molecules (including the crystallo graphic water 
molecules; 107 388 atoms): a total of 1 16 579 atoms. 

The whole system was equilibrated with velocity reassignment every 1 ps (500 steps) 
at 50 K, 100 K, 200 K and 298.15 K, for 20 ps (10 000 steps) each. Then it was equili- 
brated (without velocity reassignment) at 298.15 K for 1 ns (500 000 steps). 

The equilibration and production parts of the simulation were performed on the 
Blue Horizon, an IBM Scalable POWERparallel (SP) supercomputer, at the San Diego 
Supercomputer Center. The production part was carried out over a period of 6 months 
(not continuous). The production jobs ran on 64 to 128 processors of the Blue Horizon, 
consuming about 75 processor-months of supercomputer time. Frames were collected 
at 1 ps intervals for the production length of 5 ns, giving 5x10^ frames. 
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IILC.3 Proper radii of the passages 

The extent of opening for every passage into the gorge was characterized using a single 
variable named the proper radius, as we have introduced previously for the gorge in the 
apo-mAChE MD [2,32]. The proper radius was defined as the maximum radius of a 
spherical ligand that can go through the opening of interest from outside of the protein 
to reach the active site. This could be measured with an algorithm which finds the largest 
probe radius that generates the solvent-accessible surface with a continuous topology. 
Equivalently in practice, we detected whether the O^^ atom of Glu 202 contributes to a 
inflated solvent accessible surface that connects to the outside of the protein. The probe 
was considered capable of entry into the active site when the residues outside and those 
inside are topologically continuous. With several trials for each snapshot, we could 
narrow down the proper radius to a desired resolution. 

Three possible passages leading into the active site were located: the gorge, the 
side door and the back door. While detecting for one potential opening, regions near 
the other two potential openings were blocked to avoid contamination. The gorge was 
blocked using a sphere centered at the average position of the atoms C^, of Leu 76, 
Co of Trp 286, and Op^ of Tyr 72, with a radius of 110 pm. The blocking sphere of 
the back door was centered at the center of the 6-membered ring of Tyr 449, with a 
radius of 65 pm; that of the side door, at the Co atom of Thr 75, with a radius of 9 1 pm. 
The proper radii for the side and back doors were only calculated if their values were 
larger than 140 pm at a resolution of 10 pm. A finer resolution (5 pm) was used for 
the gorge than for the doors, and no lower bound was imposed. Therefore, a discretized 
value of the gorge proper radius p means that the true value of the proper radius p* 
satisfies p — 2.5 pm < p"" < p + 2.5 pm. For the doors, a value of proper radius p, 
when p > 140 pm, the true value of the proper radius p* satisfies p < p* < p + 10 pm 
(Figure III.3). 
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Figure III.3: Time series for the openings through the 5 ns Fas2-mAChE MD simulation. 
Top, gorge proper radius; Center, back door proper radius; Bottom, side door proper 
radius. For the latter two, points at the bottom row represent time frames with proper 
radius smaller than 140 pm. 
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IILC.4 Principal component analysis 

The C(^ coordinates of whole complex throughout the 5 ns trajectory (1 803 degrees of 
freedom: 540 C^ atoms in mAChE and 61 in Fas2) were admitted into the principal 
component analysis (PCA) [29,30]: The covariance matrix was constructed and then 
diagonalized to give the diagonal matrix of eigenvalues and the transformation matrix 
containing columns that were the eigenvectors. Using the transformation matrix, the 
projection time series along each one of the principal component was calculated. For 
each one of these principal components, the maximal and the minimal projections were 
transformed back to Cartesian coordinates. 

IILC.5 Porcupine plots 

Porcupine plots [2] may be used to show the correlation between the movement of C^ 
atoms and any desired functionally important motion. In this case, the functionally 
important motions are the extent of opening of the gorge, the side door, and the back 
door. 

The correlation coefficient between a proper radius p{t) and the x degree of free- 
dom for an C^^ atom x-(t) was defined as 

((x,.(0-(x,-),)(p(0-(p)r)), 



y((^/(0-(x,-),)')/(p(0-(p),)')^ 



(III.l) 



Similar expressions were obtained for y-{t) and ^■(?) for the y and z coordinates. These 
three correlation coefficients made up a correlation vector for the C^^ atom of each 
residue. The direction of the vector indicated where the residue was displaced when 
the proper radius became above average. 

III.D Results 

The first part of the analysis compares mAChE in the present 5 ns simulation and our 
previous 10 ns simulation [2]. The differences in structure and dynamics of mAChE are 
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investigated from rmsd plots, B factors, PCA, and average structures. The fluctuations in 
the width of the passages leading to the active site are examined and porcupine plots are 
used to pick out the motions of residues involved in their opening. In the second part 
of the analysis, the relative motion of the two proteins is assessed, and the nature of the 
interface region in the Fas2-mAChE complex is inspected by looking at the dominant 
residue contacts, particularly those that change or diverge from the crystal structure with 
time. 

III.D.l Stability of the trajectory and rmsd 

The energy components, the temperature, and the volume were inspected and found 
to have reasonable stability throughout the 5 ns simulation (data not shown), ensuring 
that the system was at equilibrium. The rmsd of the protein structures from the crystal 
structure as a function of time for both heavy atoms and C^ atoms is in Figure ni.4. The 
rmsd for the heavy atoms in the complex fluctuated but stayed around 330 pm, and for 
the Cf;^ atoms, around 280 pm. This is in contrast to the lower respective values in apo- 
mAChE of 170 pm and 120 pm. Thus, the simulated Fas2-mAChE complex appears to 
deviate more than apo-mAChE from the corresponding crystal structure. The mAChE 
part of the structure has even larger rmsd values through the 5 ns, but the fluctuation 
thereof is smaller than that of the whole complex. 

Comparing the average structure from the production run with the crystal structure 
indicates the changes that occurred during the equilibration phase. These changes were 
mainly in the termini and the segments that were flexible in the production run. The ac- 
tive site residues had also changed in its conformation during equilibration; this change 
is shown below in detail. 

IILD.2 Flexibility from B factors 

The isotropic temperature (B) factors provide a means of flnding out which parts of 
the protein are behaving differently between apo-mAChE and Fas2-mAChE. B factors 
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Figure IIL4: Root mean square deviation (rmsd) through the 5 ns Fas2-mAChE MD 
trajectory, using the crystal structure of the complex as reference. The bottom pair 
of curves are for the whole complex, with each frame superimposed over the crystal 
structure of the whole complex (IM AH). The top pair of curves are for the mAChE part 
of the complex, with each frame superimposed over the mAChE part of crystal structure 
only. In each pair of curves, the top curve is the rmsd for the heavy atoms, and the 
bottom curve is the rmsd for the C^t atoms. 
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for our MD simulations were calculated from the mean square fluctuation (msf) [14] 
(Figure IIL5). 

As in the apo-mAChE simulation, there is a wide variation in flexibility for dif- 
ferent residues. To compare apo-mAChE and the complex, the ratio of the B factors 
calculated in the complex M D trajectory over those from the apo-mAChE M D for each 
residue is calculated (Figure III.6); this information may be marked onto the protein in 
Figure IIL7. The extent of positive deviation in complex ranges from blue (least) to red 
(most). Interestingly, deviations predominate aroundanumber of residues and segments 
on the surface of the protein, namely, residues 46, 320, 462 and the segments 258-264 
(the distal small Q. loop) and 430-435. Overall, the B factors seen in the complex MD are 
larger and again support the overall trend that mAChE is more flexible in the complex 
than in the apo-form in our simulations. 

Comparing the B factors for the Fas2 calculated from this simulation and the previ- 
ous unbound Fas2 M D simulation [61], we find that the flexibility in loops I (residues 4 
to 16) and II (residues 23 to 38) that appeared in the unbound simulation has been sup- 
pressed compared to the flexibility in loop III (residues 42 to 51) andthe two turn regions 
that encompass residues 16-20 and 54-57 near the Fas2 core. 

IILD.3 Collective motion from principal component analysis 

Principal component analysis is one means of displaying collective motion. The prin- 
cipal components with the largest eigenvalues represent the largest scale motions. The 
calculated eigenvalues decrease in magnitude quickly but smoothly: the 8th largest is 
less than 0.1 of the largest, and the 38th largest is less than 0.01 of the largest. The 
Cartesian-coordinate structures for the maximal and minimal projections of the princi- 
pal components with the two largest eigenvalues (1.58 nm^ and 0.87 nm^) are shown in 
Figure IIL8. 
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Figure 111.5: Top, B factor for each residue in the mAChE part of the complex; bottom, 
that in the Fas2 part. Dashed line, B factor of the Cc^ atoms as reported in the crystal 
structure (1 KU6); solid line, calculated for the whole residue from the msf in the 5 ns 
Fas2-mAChE M D trajectory. 
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Figure IIL6: Ratio of B factors for each residue in mAChE, calculated as that from 
Fas2-mAChE MD over that from the apo-mAChE M D . The dashed line marks 1. 



IILD.4 Comparing the average structures of mAChE 

For both the apo-form and the complex simulations, we produced the average structure 
of mAChE by superimposing the enzyme conformations onto that in the crystal struc- 
ture and taking the mean of the atom positions. The resulting two average structures 
are displayed in Figure IIL9. The distances between the C^ atoms, and the absolute 
values of the <p and I//" angle differences for every residue in each average structure were 
calculated (Figure III. 10). When comparing these average structures, we have to keep in 
mind that many differences in the distances and angles might be linked to the differences 
in crystal packing and in the resolutions of the starting crystal structures. 

While the overall agreement between the two structures is close, there are many 
places where they differ in small ways, and in a few loops the difference is quite sub- 
stantial, such as residues 258 to 264 (the distal small Q. loop), 370 to 393, and 436 
to 442, as well as the C-terminus. A greater variation is seen in the dihedral angles, but 
evidently many of these changes appear to cancel out and lead to little net movement in 
Ca position. 

The active site arrangement was changed in the Fas2-bound complex, with critical 
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Figure III. 7: Ratio of B factor for each residue in mAChE, calculated as that from Fas2- 
mAChE MD over that from the apo-mAChE M D , colored on the mAChE structure; each 
sphere represents one residue. The ratio decreases from red to blue. The active site 
Ser 203 is shown in the space-filling representation. Generated using OpenDX (IBM, 
White Plains, New York) with Chemistry Modules [60]. 
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Figure 111.8: Structures representing the extremal projections on the principal compo- 
nents with the largest eigenvalue, 1.58 nm^ (left) and the second largest eigenvalue, 
0.87 nm^ (right). Red, maximal projection; blue, minimal projection. The thicker back- 
bone is the Fas2 part of the structure. Generated using VMD [62]. 



implications for the function of mAChE. The change was seen by comparing the average 
structures from the apo-form and the complex trajectories (Figure III. 11): In the apo- 
form average structure, the arrangement was conducive to proton transfer placing the 
imidazole ring of His 447 between the carboxyl group of Glu 334 and the hydroxyl 
group of Ser 203. Visual inspection of the complex trajectory and the average structure 
therefrom showed that by many dihedral angle changes in the backbone and side chains, 
the His 447 imidazole ring moved to a conformation where it is almost orthogonal to 
the route of proton transfer. In addition, the side chain carboxyl group of Glu 334 has 
moved away from the route. This disruptive change in conformation occurred in the 1 ns 
equilibration procedure before the 5 ns production run, and persisted through the 5 ns 
production trajectory. 
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Figure III. 9: Comparing the average structures from the 10 ns apo-mAChE (red) and 
from the 5 ns Fas2-mAChE (green) trajectories. The backbones are displayed in the 
ribbon representation, with the a-helices and j8-sheets shown. Residue Ser203 from the 
average structure of the apo-mAChE trajectory is shown in the space-filling representa- 
tion. Generated using Insight 11. 
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Figure III. 10: Comparing the average structures from the 5 ns Fas2-mAChE and the 
10 ns apo-mAChE trajectories. Top: the distances between the C^^ atoms of each residue 
in the two average structures. Center: the absolute values of the 0-angle differences. 
Bottom: the absolute values of the l//-angle differences. 
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Figure III.l 1: Catalytic triad of mAChE. The average structures from both simulations 
are shown, after all-atom superimposition of the mAChE molecule. Gray, the apo- 
mAChE MD; colored, from the Fas2-mAChE MD. From left to right, the residues are 
Ser 203, His 447, and Glu 334. Generated using Insight 11. 

IILD.5 Opening of the passages 

The one entrance into the active site that has been observed in all available crystal struc- 
tures is the main active site gorge itself. The distribution for the gorge proper radius in 
the 5 ns Fas2-mAChE simulation is shown in Figure III. 12; the time series thereof is in 
Figure III. 3. As in the 10 ns apo-mAChE MD, we found two distinct local maxima in 
the probability distribution of the proper radius [32]. This affirms that there are indeed 
two states in mAChE gorge fluctuation: one narrower and one wider. In comparison, 
with Fas2 bound, the narrower state becomes more favored than in the apo-form simula- 
tion, and the distribution as a whole moves towards a smaller value for the gorge proper 
radius. 

In addition, we produced the selective average structure for only the frames with 
an open gorge (proper radius greater than 127 pm) and that for only the frames with a 
closed gorge. The selective average structures for the open and closed version do not 
appear to be significantly different. This is also true for the selective structures of the 
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Figure 111.12: Distribution of the gorge proper radius. Solid line, that in the 10 ns 
apo-mAChE trajectory; dashed line, twice that in the 5 ns Fas2-mAChE trajectory (so 
normalized for ease of comparison with the 10 ns data). 

side door, and those for the back door. (For these two, an opening is defined as having 
a proper radius greater than 140 pm.) If we focus on the gorge, the distance between 
the Ca atoms of Tyr 124 and Phe 338 that straddle the gorge in the selective average 
structure with an open gorge is 1.288 nm; the distance in the structure with a closed 
gorge is 1.281 nm; the difference is only 7 pm. For in the apo-mAChE trajectory, the 
same distance measurement has 1.485 nm in the open gorge case, and 1.440 nm in the 
closed case; the difference is 45 pm and is comparable to the difference in gorge radius 
between the two states seen in Figure 111.12. 

A much larger and possibly functionally important difference was observed for 
the back and side doors. In the apo-form simulation, the back door was open to water 
molecules less than one hundredth (0.01) of the frames and the side door was never 
observed to be open. However, in the complex both side and back doors are observed 
to be open for a significant amount of simulation time. As shown in Figure III.3, 917 
out of the 5 000 frames collected (0.18) have a back door opening event (proper radius 
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larger than 140 pm); 664 out of the 5 000 frames (0.13) have a side door opening event 
(defined likewise). There are 1 158 frames (0.23) with either of the alternative passages 
open, and 423 frames (0.08) with both open. 

To probe the access to the active site by ligands or solvent molecules through the 
gorge when Fas2 is bound, we blocked the side and back doors, and measured the largest 
opening in the gorge region with Fas2 in place. The gorge is open to a spherical probe 
of radius 140 pm in 233 frames (0.05), while in only 9 frames is the gorge open for a 
spherical probe of radius 160 pm. 

IILD.6 Motions of passages using porcupine plots 

Porcupine plots are an alternative means of displaying concerted protein motion. In the 
present work, they are used to extract out protein motion correlated with the radius of 
each passage found into the gorge. Three porcupine plots, for the gorge, the side door, 
and the back door, are shown in Figure III.l 3. The coloring indicates the extent to which 
a residue's motion is correlated with the opening of the given passage, with red the most 
and blue the least. 

Different parts of the protein are seen to be associated with each passage. Among 
the largest vectors in the side door porcupine plot is at residue 86, whose C^ atom is 
close to the side door region. Residues 445 and 446, in the back door region, have two 
of the largest vectors in the back door porcupine plot. 

III.D.7 Movement and contacts of Fas2 

The first issue addressed here is the extent of Fas2 motion relative to mAChE. The 
trajectory of the center of geometry of Fas2 (defined as the average position of the C^ 
atoms) relative to that of mAChE was calculated. It was found that this Fas2 trajectory 
is contained within a sphere of radius 250 pm. Thus it appears that the two proteins 
remain largely fixed with respect to each other for the duration of the simulation. 

The second question is what residue contacts are made at the interface and whether 
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Figure 111.13: Porcupine plots. A, for the gorge in the apo-mAChE trajectory; B, for the 
gorge in the Fas2-mAChE trajectory. In these two, the same coloring scheme is used. 
C, for the side door in the Fas2-mAChE trajectory; D, for the back door in the Fas2- 
mAChE trajectory. In each, three residues are shown in the space-filling model: Ser 203 
(marking the bottom of the gorge and the active site), top; Tyr 449 (marking the back 
door), right; Leu 76 (marking the side door), bottom left. Generated using OpenDX with 
Chemistry Modules [60]. 
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any of these change. To pick out interesting contact pairs at the Fas2-mAChE complex 
interface, we record any distance between a heavy atom in Fas2 and another in mAChE 
that is shorter than 500 pm formore than half of the 5 x 10^ frames in the 5 ns trajectory. 
There are 599 contact pairs between heavy atoms in the interface. Eighty-three (83) of 
these are polar (between non-carbon heavy atoms) contact pairs; on the residue level, 
there are 39 interesting polar residue-residue contact pairs. For the 83 polar atom- 
atom pairs, we tallied the distributions of their distances through the 5 ns trajectory. If 
there is more than one maximum to a 50 pm resolution in a distribution, the residue 
can be considered more dynamic; otherwise, it is considered static, as the distance only 
fluctuates around one value. The details of these pairs are shown in Figure III. 14 and in 
the Supporting Information of [3]. 

III.E Discussion 

Fas2 may influence the function of mAChE in a variety of complementary ways: Fas2 
may sterically obstruct the entrance of the mAChE gorge, allosterically change the con- 
formation of parts of the enzyme, affect the fluctuation behaviors of the openings dy- 
namically, or affect the efficiency of the catalysis at the active site, among other things. 
Hence our results will be discussed in light of these effects that may conceivably be 
adopted in natural design of Fas2. 

III.E.l Steric obstruction of the gorge entrance 

As the eigenvalues from the PCA decrease rapidly in magnitude, the largest collective 
motions are represented in the few principal components with the largest eigenvalues. 
Looking at two of these components (Figure III. 8), we do not see any collective move- 
ment that suggests Fas2 moving away significantly from the binding configuration where 
it sterically obstructs the gorge entrance. The Fas2 center of geometry movement rel- 
ative to mAChE was also small; indeed, this motion can be contained by a sphere of 
radius 250 pm. 
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Figure III. 14: Interface contact pairs. Top, structure of the residues participating in 
contact at the start of the 5 ns production trajectory (at ns; that is, right after the 
1 ns equilibration); bottom, that at 5 ns. The residues shown are those participating in 
interfacial interactions more often than not during the 5 ns, both polar and nonpolar. 
Residues belonging to Fas2 are shown in green; those of mAChE, blue. The yellow 
bars represent polar interactions with only one maximum in the distance distribution 
(static); the red bars represent polar interactions with more than one maximum (dy- 
namic). Generated using OpenDX with Chemistry Modules [60]. Inset, the backbone of 
the relevant segments in Fas2 (gray) and in mAChE (white), with key residues labeled 
(Fas2 in parentheses); generated using Molscript [56]. 
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With Fas2 bound, the active site is rarely accessible through the gorge to allow 
entry of a spherical ligand with radius 160 pm (much smaller than the acetylcholine 
profile). Since it is nearly impossible for a ligand this size to move through the gorge, 
access to the active site is only likely through the alternative passages; in this simulation, 
we observed a large increase in their opening likelihood compared to the apo-mAChE 
trajectory. 

IILE.2 Interactions at the interface 

The polar interface contact pairs are listed in the Supporting Information of [3]. The par- 
ticipating mAChE residues mostly fall into two groups: the Q. loop group (residues 69 
to 96) and the residues near the peripheral site (roughly from Asp 283 to Ser 293, and 
Tyr 341). On the Fas2 side, the participating residues are mostly from loop I or II; only 
Asn 47 is from loop III, and Tyr 61 is at the C-terminus of Fas2. It is notable that nearly 
all Fas2 residues at this interface participate in polar interactions (Figure III.14). 

The contact regions of Fas2, namely loops I and II, have higher or comparable 
flexibility compared to loop III in the previous unbound Fas2 simulation [61]. When 
bound to mAChE, they are surpassed in flexibility (measured by the B factors) by the 
non- contacting loop III. 

Except for the few contacts made by the Fas2 core region, we do see most of the 
contacts reported in the crystal structure [28], and when we do not see the exact inter- 
action, usually a neighboring residue has a contact pair. In short, loops I and II of Fas2 
have contacts with the mAChE Q. loop, while only loop II is in charge of interacting 
with the peripheral site. The flexibility of these contact loops is decreased relative to 
other parts of Fas2 upon binding. 

IILE.3 Allostery: mAChE conformation and flexibility changed 

As we embark in comparing the apo-mAChE and the complex trajectories, there is a 
caveat we have to point out first: With nanosecond-scale simulations, the trajectories do 
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not sample adequately all conformations available for a protein of the size of mAChE 
to explore. Some of the differences observed may be real and indicative of changes 
in protein functions under different circumstances; others may simply be trajectory- 
specific idiosyncracies that are imprudent to generalize. 

Most residues in mAChE exhibited larger B factors in the Fas2-mAChE simula- 
tion than in the apo-mAChE one (Figure III.6). In particular, the distal small Q. loop 
from Pro 258 to Gly 264 had an extremely large increase in B factor as compared to 
that in the apo-mAChE simulation. Notably, this loop could not be resolved in the 
Fas2-mAChE complex crystal structures (IMAH and lKU6) [28] because of its high 
flexibility, while in the tetrameric mAChE structure (IM A A) [57], where it participates 
in the crystalline tetramer interface and has stabilizing interactions, it was resolved in 
two of the monomers. The flexibility of this loop, as observed in the simulation, may be 
the reason for its low resolution in the complex structure. 

Our previous apo-mAChE simulation [2] used the mAChE molecule from the Fas2- 
mAChE complex structure at 0.32 nm (3.2 A) resolution (IMAA); the present Fas2- 
mAChE simulation used a Fas2-mAChE structure at 0.25 nm (2.5 A) resolution, a res- 
olution improvement that would be expected to generate lower B factors for the Fas2- 
mAChE simulation. Yet, we observed a higher flexibility in the complex trajectory than 
in the apo-mAChE trajectory. 

By visual inspection, the most impressive differences between the mAChE struc- 
tures from the apo-mAChE and Fas2-mAChE complex trajectories are in the positions 
and flexibilities of the distal small Q. loop from Pro 258 to Gly 264 (at the leftmost of 
Figure III.9) and the remote segment around residue 380 (at the bottommost of the same 
figure). Despite being far from the Fas2 binding site, the distal small Q. loop seems 
to acquire a different average position and become more flexible (as seen in the B fac- 
tors) upon Fas2 binding. This is in stark contrast with fluorescence spectroscopy results 
where the fluorophore attached at residue 262 of this distal loop indicated little alter- 
ation in flexibility [63]. To reconcile this discrepancy, we recall that the distal loop was 
built in artificially in the two M D simulations. Thus it may not be justified to compare 
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the difference in flexibility in MD of this loop, considering its local lack of resolution in 
the crystal structure compounded with the ambiguity introduced by the artificial model- 
building process. 

IILE.4 Dynamical inhibition of gorge opening 

Contrary to previous expectations [28,57], the observed change in the fluctuation (as 
indicated by the B factors) and conformation of the Fas2-bound mAChE Q. loop are 
only moderate and not outstanding in view of other parts of the molecule. Evidences 
both from simulations and from fluorescence spectroscopy experiments [55] indicate 
that, though mAChE and Candida rugosa lipase are related in their a/j8 hydrolase fold, 
the former does not seem to have its open and closed states characterized by the well- 
defined hinge motion of the Q. loop as the latter [64]. 

Considering that the two peaks in the probability distribution function of the gorge 
proper radius were separated by about 50 pm in the apo-mAChE trajectory [32], it is 
not surprising that the difference between the Tyr 124-Phe 338 C^^ distance is 45 pm. 
On the other hand, the smaller 7 pm difference in the complex trajectory manifests 
that the bottleneck region of the gorge is no longer coupled to the Tyr 124-Phe 338 
distance upon Fas2 binding, and implies significant change in the mechanism of the 
gorge opening behavior. 

The binding of Fas2 indeed favors a narrower gorge size, as shown in Figure III. 1 2. 
Not only does the whole distribution shift towards a smaller value of the proper radius, 
but the narrower of the two distinct states becomes more likely to occur in the complexed 
trajectory. Comparing the porcupine plot for the gorge with that from the apo-mAChE 
MD [2,65], we find that the concert of motions near the Fas2 interface seems to be 
repressed in the Fas2-mAChE trajectory. Such effect of dynamical inhibition [26] is 
impressive, though the mechanism by which Fas2 imposes this effect is still not clear to 
us. 
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IILE.5 Changes in the alternative passages to the active site 

In contrast to the apo-mAChE trajectory where only one hundredth (0.01) of all the 
frames collected gave an open back door [2], we have an increase in the likelihood of 
opening to 0.18 when Fas2 is bound. In addition, the side door also has a 0.13 chance 
of being open in the complex trajectory (Figure IE. 3). In the porcupine plots, most 
concerted motions show up around the respective alternative passages (the side door and 
the back door). These increases in the number and likelihood of alternative openings, 
as previously suggested, may indeed be responsible for some of the residual, or even 
enhanced, substrate catalysis in Fas2-bound AChE [28, 50-52]. 

In a previous MD simulation of mAChE complexed with the active site ligand hu- 
perzine A [66], opening of the side door is observed more frequently than in the apo- 
form simulation [27]. In another simulation, where the T. californica AChE dimer is 
complexed with the active site inhibitor tacrine [25], side door opening is also observed. 
It is in this context that we report here the significant increase of side door opening 
events induced by an inhibitor that does not bind at the AChE active site. This im- 
plies that alternative passage opening events may be promoted by either active site or 
peripheral site binding. 

IILE.6 Changes in the conformation of the active site 

Comparing the average structures (Figure III.l 1) from the apo-mAChE and Fas2-mAChE 
trajectories, we notice dramatic differences in the mAChE catalytic triad: Upon Fas2 
binding, the functional groups in Glu 334 and His 447 have moved out of the strate- 
gic positions conducive to proton transfer that are observed in the apo-mAChE M D and 
the crystal structure. Such movements are achieved through movements of the back- 
bone and several dihedral angle changes in the side chains, and do not appear to be 
trivial. Specifically, in the average structure from the Fas2-mAChE MD, the Glu 334 
carboxylic group now points away from the neighboring His 447 in the complex M D , 
and the His 447 imidazole ring is now at an angle from its position in the apo-mAChE 
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MD, out of alignment for proton transfer. 

Crystallography of T. californica AChE complexed with the organophosphorus in- 
hibitor VX showed a movement of the active site histidine [67]. Monoclonal antibodies 
that bind to the peripheral site has been reported to allosterically affect the orientation 
of Trp 86 near the active site [68]. In addition, kinetic data have suggested that Fas2 
inhibits AChE by disrupting the conformation of the active site, thus slowing down the 
proton transfer steps [49]. Our observation here is consistent with such a mechanism. 
However, as the Fas2 binding site does not neighbor the active site, the sequence of 
mechanical processes through which Fas2 exerts its influence is not immediately osten- 
sible. Even though more rigorous procedures (such as the use of particle-mesh Ewald 
for the long-range electrostatics) were used here than the early M D work on T. califor- 
nica AChE where active site disruption was relieved by addition of absent counterions in 
the gorge [69], the changes in the active site presently observed can only be considered 
suggestive at this stage. 

III.F Conclusions 

In the Fas2-mAChE complex, Fas2 binds to the mAChE Q. loop and peripheral site with 
excellent surface complementarity using many polar and hydrophobic interactions from 
loops I and II [28]. Our 5 ns simulation of this complex shows that Fas2 binding dynam- 
ically disfavors the opening of the mAChE active site gorge but increases the likelihood 
for back door and side door opening, a feature that may contribute to the residual (and 
occasionally enhanced) catalytic activity of the Fas2-bound enzyme observed in solu- 
tion. As in the apo-mAChE simulation, the fluctuating behavior of the mAChE gorge 
in the Fas2-bound complex continues to present complexity, though selective averaging 
suggests a change in the location of the bottleneck region. Fas2 binding also increases 
the flexibility of two surface segments remote to the binding site, and, most impres- 
sively, disrupts the catalytic triad arrangement of the mAChE active site. In sum, we 
have observed steric, allosteric, and dynamic effects likely to represent the components 
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of the mechanism through which Fas2 influences the function of mAChE. 
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IV 

Finite element simulations of 
acetylcholine diffusion in 
neuromuscular junctions 

What I cannot create, I do not understand. [70] 

IV.A Abstract 

A robust infrastructure for solving time-dependent diffusion using the finite element 
package FEtk has been developed to simulate synaptic transmission in a neuromuscular 
junction with realistic postsynaptic folds. Simplified rectilinear synapse models serve as 
benchmarks in initial numerical studies of how variations in geometry and kinetics relate 
to endplate currents associated with fast-twitch, slow-twitch, and dystrophic muscles. 
The flexibility and scalability of FEtk affords increasingly realistic and complex models 
that can be formed in concert with expanding experimental understanding from electron 
microscopy. Ultimately, such models may provide useful insight on the functional im- 
plications of controlled changes in processes, suggesting therapies for neuromuscular 



This chapter is reproduced with permission from the manuscript Tai, K.; Bond, S. D.; MacMillan, H. R.; 
Baker, N. A.; Hoist, M. J.; McCammon, J. A. submitted to Biophys. J. 2002. 
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diseases. 

IV.B Introduction 

The neuromuscular junction (NMJ) is the point of communication between neurons 
and muscle fiber in the orchestration of muscle contraction. This study encompasses 
the release of neurotransmitter acetylcholine (ACh), its hydrolysis with acetylcholin- 
esterase (AChE) clusters, and its reactive presence near acetylcholine receptor (AChR) 
molecules. Past computational modeling of synaptic transmission on this scale has 
involved either differential equations governing continuum reaction-diffusion [71,72] 
or particle methods such as Brownian dynamics and Monte Carlo [73,74]. Here, we 
present an improved finite element infrastructure to solve continuum reaction-diffusion 
using FEtk [75], an efficient platform for adaptive multiscale modeling. These efforts 
highlight the flexibility of finite elements to evolve with improved understandings of re- 
action kinetics. Now, with the capacity to represent realistic NM Js, comparative studies 
can be conducted with the guidance of coordinated experimental data. Ultimately, this 
could provide insight on the functional implications of NMJ variations associated with 
neuromuscular diseases, such as muscular dystrophy and myasthenia gravis. 

IV. C The neuromuscular junction 

A typical NMJ features a smooth presynaptic neuron membrane and a folded postsynap- 
tic muscle surface of crests and troughs. When an action potential reaches the end of the 
nerve, it results in a localized influx of calcium ions that, in turn, causes vesicles to fuse 
to the neuron membrane and release ACh. In our model, vesicles are treated as having 
just opened and we do not yet account for subsequent relaxation of the neuron mem- 
brane. Eventually, this effect may be incorporated into FEtk, along with spatial control 
of vesicle placement according to the varying presence of calcium ions. 

ACh diffuses across the synaptic cleft and potentially binds to acetylcholine recep- 
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tors, ion channels embedded in the postsynaptic folds. AChR ion channels are multi- 
protein membrane- spanning complexes found at packing densities of up to 10 000 fim~^ 
at the crests of the postsynaptic folds. In the absence of ACh, an AChR channel is imper- 
meable to ion flow. However, once two ACh molecules bind to it, AChR opens and ions 
flow through the muscle cell membrane: sodium inward, potassium outward. This ion 
flow defines an endplate current (EPC) across the muscle membrane that, when strong 
enough, induces contraction. Roughly 1 ms after activation, ACh molecules are released 
back into solution and AChR ion channels close. 

ACh molecules remain in the cleft until they are hydrolyzed by AChE, the biomolec- 
ular "off-switch" for synaptic transmission. AChE is present as clusters of three tetramers 
suspended by collagen stalks bound to the muscle membrane at varying density (from 
600 jUm~^ to 2500 /xm"^) throughout the postsynaptic folds. As an extremely fast 
enzyme capable of destroying ACh molecules at rates approaching theoretical limits, 
AChE provides a very efficient mechanism to terminate synaptic transmission for sub- 
sequent signaling [1,9]. 

Experimental data are available on the ultrastructure and activity of NMJs of dif- 
ferent muscle types to guide the initial development of mathematical models [76-81]. 
Ultrastructural differences between the NMJs in vertebrate fast (twitch or extensor digi- 
torum longus) and slow (tonic or soleus) muscles have been observed [10, 82-84]. Also, 
geometric and reactive deviations in mouse NMJs due to muscular dystrophy have been 
documented [11, 85-87]. Local measurements of miniature endplate current (mEPC) are 
used to infer the functional implications of such structural differences on an NMJ's ef- 
ficiency. Having established our simulation infrastructure, we can begin to recreate the 
effects of the various parameters of different muscle types in silica. It is evident that 
more experimental data coordinated with simulations are needed to bring increasingly 
realistic models to maturity. 
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IV.D Mathematical setting 

IV.D.l Continuum formulation 

In a continuum model, the concentration of ACh, u{x,t), is assumed to satisfy the time- 
dependent diffusion equation in the synaptic cleft, Q., with appropriate conditions on its 
boundary, dQ.. Formally, the concentration varies from its initial state, m(x,0), according 
to 

^ ' ^ -V-DVu{x,t) = in^, (IV.l) 



6t 



n(x)-DVu(x,t) = <^ , (IV.2) 

on do. -do.. 



-act 



where D = 4.0 x 10"'* firn^ -/is"' is the diffusion coefficient, n(x) is the outward unit 
normal, and f^^^ct ^ ^^ denotes the cumulative reactive surface of AChE with specific 
reactivity k^^^ [71]. Note that Equation IV.2 states a zero-flux condition on all bound- 
aries except the surfaces representing AChE clusters, dO.^^.^, for which a linear reaction 
scheme yielding a Robin boundary condition is now assumed. Estimation of the specific 
reactivity, k^^^, is described in each example. Attempts to lift this linear assumption on 
AChE binding and use more involved boundary conditions will be the focus of future 
research. 

IV.D.2 Discrete formulation 

We use the backward Euler algorithm to define a method of lines and reduce Equa- 
tions IV.l and IV.2 to the following elliptic problem at time /„, knowing the previous 
concentration, u{x,t^_^ ): 

u(x,tn) ~ m(x,C I ) 
-V-DVu{x,tn) + ^^^ ^ "~^^ =0 inn, (IV.3) 

n{x)-DWu = { ''' ''' . (IV.4) 

ondQ.-dQ.^,i 
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It is convenient to write w„ := uix^t^ ) and define 

b{x,un,u^_,) := ^ '" _ ^ '"-" . (IV.5) 

Integrating by parts, we obtain the weak form of Equation IV.4, 

/ £)Vw,j-Vvdx— / vDVun-ndS+ / /?(x,w,j,w , )vdx = 0, 
Jo. Jsa Ja 

Vv e V, (IV.6) 

where V is the test space [21, 75]. Enforcing boundary condition on dQ., the weak form 
becomes 

/ Z)Vw„ • Vvdx + ^' . / UnvdS+ I ^(x,w„,M _,)vdx = 0, 
Jo. -^da^ci ^^ 

Vv e V. (IV.7) 

For a discrete solution to Equation IV.7, a finite element space V = span{(^| , . . . , <^^} C 



V is used. The Galerkin approximation 



A^ 



(=1 
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(IV.8) 



satisfies the discrete weak form 



j DVul-V^.dx + k[,J «i>,.d5+/fo(x,«^,«tl)<^,dx = 0, 

V(^,.e{(^i,...,(^^}, (IV.9) 

knowing the discrete solution from the previous timestep, wjj_^ = Y.i=\ ^°i^i ^ ^ • To 



formulate Equation IV.9 into a matrix equation, we write the terms 



It follows that 



f ^ \ f 



n-l / ( = 1 L 



Au+— l\/l(u-u°) + /c' .Fu = 0, 

At 



(IV. 10) 
(IVll) 
(IV. 12) 

(IV. 13) 



where the stiffness matrix A = [A. ] = J^DV^- • V^- dx , the mass matrix M = [M--] 
Iq. ^i^i dx , F = [i^. -] = /^^ ^1^1 dS , and the solution vectors u = [wj, u° = [uf]. 
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IV.D.3 Numerical solution and visualization 

Several methods of lines to solve Equations IV. 1 and IV.2 have been built in the context 
of FEtk [75], a finite element package written in "clean object-oriented" C language. De- 
tails on the resulting linear systems at each timestep, using the backward Euler method 
of lines described above. At each timestep, the conjugate gradient method is used with 
termination chosen such that time-truncation error is less than 10"^*^. Given a surface 
triangulation of an NM J boundary with realistic folds, a robust and flexible mesh gen- 
eration package, NETGEN [88], develops a structured volume mesh that has been inter- 
faced with FEtk (Figures IV.l and IV.8). A structured mesh fine enough to sufficiently 
describe the constraining surface geometry is used in the examples that follow, and re- 
fining this mesh (i.e., halving edge length) has no appreciable effect. Similarly, the 
timestep has been chosen to adequately resolve a postsynaptic response curve; e.g., the 
range 10~' fis to lO' fis was studied to ensure that the simulation is not overly sensitive 
to the chosen timestep, 1 fis. Computing 1 000 timesteps for a system of approximately 
33 000 vertices (such as the rectilinear synapse model) takes less than 1 hour on an 
Intel (Santa Clara, California) Xeon 1.80 GHz personal computer. The classical duo 
of backward Euler method of lines combined with conjugate gradient is nearly optimal 
for ordinary diffusion. Moreover, FEtk memory usage scales linearly with the number 
of vertices. The FEtk output formats are readable using MATLAB (MathWorks, Natick, 
Massachusetts), GMV (Los Alamos National Laboratory, New Mexico) and OpenDX 
(International Business Machines, White Plains, New York). 

IV.D.4 Postsynaptic detection and AChR binding 

Equations IV.l and IV.2 quantify the presence of ACh in the NM J, but relating this to 
experimentally observable mEPC involves estimating not only the density of AChR, but 
also the stochastic nature of its binding and retention. By defining the postsynaptic 
detection level, L(t), as the weighted surface integration 

L(0 := / Y^{x)u(x,t)dS, (IV.14) 

■J do. 
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where 7j^(x) is AChR density, we only observe general timescale trends of observed 
mEPC, such as rise time and signal duration. 

Toward more accurate estimations of mEPC , we account for ACh-AChR binding in 
the model to properly estimate the time course of open AChR ion channels. This entails 
appealing to the equilibrium mechanism 2ACh + AChR^^i^^ggj ^ (ACh)2(AChRopen). and 
assuming that AChR is saturated with ACh. For the equilibrium constant, we have 
K = y^ 6 /[u^y^{\ — 0)], where 6 is the fraction of open AChR channels. This implies 
Q = Ku^ / (1 + Ku~ ), leading to a new measure of postsynaptic response 

A(t) := / 7r (x) I \, dS. (IV.15) 

According to [89], K = 3.6x 10" mM~^ = 1.0 x 10~^ fim^. While A(?) is an improve- 
ment, it still does not consider the effect of ACh retention when bound to AChR. To 
do this is to alter the zero-flux conditon on the postsynaptic folds with an integral that 
records accumulations of ACh bound to the ion channels. This sophistication will soon 
be incorporated into our modeling framework for realistic NM Js so that proper bench- 
marking with coordinated mEPC observations can be conducted. For brevity, we do not 
present plots of A(?) in all the examples that follow, as they cannot yet be compared 
with experimental data. 

IV.E Neuromuscular junction models 
IV.E.l Rectilinear benchmarking 

As a first point of contact with existing literature, we have created a simplified recti- 
linear NMJ model similar to that in [71]. We introduce this model before presenting 
preliminary studies on the functional implications of variations in this geometry. All 
simulations are in micrometers and microseconds, with concentrations in particles per 
volume and densities in particles per area. Rectilnear models have also been studied 
withMCell[73,74]. 
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Neuromuscular geometry and initial conditions 

In [71] a simplified rectilinear model of an NMJ is considered with 3 identical secondary 
clefts and the following measurements (Figure IV.l): primary cleft length = 2.0 /xm; 
primary cleft width = secondary cleft width = 0.05 fim; separation between secondary 
clefts = 0.5 jUm; secondary cleft depth = 0.8 fim. 

We include a vesicle as a sphere of radius 0.024 fim fused to the middle of the 
presynaptic membrane. Its center is placed 0.016 fim above this membrane, leaving a 
circular area of radius 0.018 fim as the pore opening. The initial ACh concentration is 
only nonzero inside the vesicles, where it is 300 mM = 1.8 x 10^ fim~^ [5,71,90,91]. 
With the vesicle size stated above, this translates to 6.06 x 10^ ACh molecules present. 
We recognize that this is a rather important yet simplistic aspect of our simulation, and 
sensitivity to the initial distribution of ACh will be included in future studies. 

AChE clusters are treated as three-dimensional "holes" in the NMJ mesh with 
boundary condition Equation IV.2. The previous finite element simulation [71] uses an 
AChE cluster density of about 240 jUm~^ . By comparison, in [73], an active site density 
of 1800 jUm~^ in the primary cleft (equivalent to a cluster density of about 150 fim~^) 
and twice that in the secondary cleft is used. The length of the collagen stalk that con- 
nects AChE from the synaptic basal lamina is on the scale of 0.05 /xm, depending on 
the species [92]. As done previously, we place AChE midway between the presynaptic 
and postsynaptic membranes in our rectilinear models. We have chosen to use cubic 
boxes with 0.02 fim sides spaced 0.1 jUm apart in a square array spanning the primary 
and secondary clefts to afford a cluster density of 100 fim~^ . 

Reactivity of AChE 

The reactivity of an AChE cluster in our model is relative to the surface area used to 
represent it. Following [71] we determine k[^.^ by conducting several separate but simple 
simulations. We represent an AChE cluster as a single cube with 0.02 fim edges, and 
place it at the center of a cube with 0.30 fim edges. The boundary condition for ^^^ct 
in Equation IV.2 is used on the surface of the interior cube, the AChE cluster. On the 
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Figure IV. 1 : Three views (edges, outside, inside) of the finite element mesh for the recti- 
linear model of the neuromuscular junction, with three secondary clefts and a spherical 
vesicle fused to the presynaptic membrane. The cubes represent acetylcholinesterase. 
1°, primary cleft; 2°, secondary cleft. 
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exterior cube boundary dQ. — dQ.^^^, the concentration of ACh is held at a different 
constant value for each "trial" simulation. The "current" of ACh consumed by AChE is 
defined to be 

J := [ h-DVudS. (IV. 16) 

Naturally, / must also equal the flux of ACh into the domain due to imposing 
constant ACh concentration on the outer boundary. Assuming a nonsaturated steady 
state and the catalytic scheme E ^ S ^^ E ^ P, the flux into the domain is estimated 
in linear proportion to the concentraion of ACh maintained on the boundary [71,93]. 
Writing this proportionality constant as k, we set 

and consider simulations with different sustained steady state levels of ACh on the outer 
boundary. Since k = k^^jKy^ ^ 3 x 10^ M~' -min"^ = 5 x 10~^ jirn^ • iis~^ has been 
determined by experiment [52], enforcing Equation IV. 17 for each steady-state leads 
to an estimation of k'^^^. FEtk with the backward Euler timestepper is used to compute 
steady state ACh diffusion between the concentric cubes, with validation from an ex- 
plicit steady-state solver. Sampling the range 1.00 x 10^ jUm~^ to 1.00 x 10^ jUm~^ 
(that is, 0.166 mM to 16.6 mM) as the imposed ACh concentration on the outer bound- 
ary, we find k[^.^ = 2.0 x 10~^ fim • fis~^ . This is slightly larger than the value used 
previously [71], due to the fact that we are using here a different area for the active site 
of AChE. With the enzymatic activity of AChE approximated. Equations IV. 1 and IV.2 
can now be solved for the rectilinear N M J . 

Solution and postsynaptic detection 

The time course of ACh diffusion in the rectilinear NM J model has been solved using 
the backward Euler and FEtk (as outlined above) with timesteps of 1 fis. Figure IV.2 
shows the total number of ACh molecules over 400 jUs; the same decay has been shown 
in the solid curve in Figure 3 of [71]. 
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Figure IV.2: Total number of ACh molecules in the rectilinear NMJ model decreases 
over time. 

Determining the postsynaptic detection level entails defining AChR density. MCell 
models typically use a graduated receptor density function 7j^ [73,74]. From the crest 
to 0.2 fim to 0.3 jUm into the fold, y^ ranges from 7 000 fim~^ to 10 000 fim~^ . 
It is 2 000 fim~^ to 3 000 fim~^ another 0.2 fim to 0.3 fim from the crests, be- 
low which density falls off dramatically. Here we use a piecewise constant density of 
Y^ (x) = 8 500 jUm~^ from the crests down 0.25 fim into the secondary clefts, y^ (x) = 
2 500 jUm~^ an additional 0.25 jim into the secondary clefts, and then zero down to 
the troughs. The resulting postsynaptic detection level L{t) is shown in Figure IV.3. 
The "rise time" (duration from the release of the vesicle to the peak) is about 20 fis, 
much shorter than that of mE P C . However, since we do not yet consider ACh binding of 
AChR, this is reasonable. Moreover, it is comparable to the previous observation [71]. 

IV.E.2 Variations in muscle type 

Fast-twitch and slow-twitch neuromuscular junctions 

Simple rectilinear NMJ models for fast-twitch and slow-twitch muscle have been built 
according to the primitive dimensions listed in Table IV. 1. Each model contains one 
vesicle, three identical secondary clefts, and the AChR distribution stated above. Also, 
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Figure IV.3: The postsynaptic detection level L(t) in the rectilinear NMJ model. 

normal dystrophic 

fast slow fast 

primary cleft width / fim OJO OlO OlO 

secondary cleft depth / fim 1 .00 0.75 0.50 

secondary cleft separation / fim 0.25 0.50 0.75 

secondary cleft width / fim 0.10 0.20 0.10 

Table IV. 1: Geometric dimensions for the recilinear models of different muscle types 
[10,83,86]. 



as before, AChE clusters are 0.02 fim edged cubes raised 0.025 jim above the postsy- 
naptic membrane with density 100 /im~^. 

In Figure IV.4, the rise time of L(?) for the fast-twitch NMJ model is 45 /is, and a 
peak value of 4.25 x 10^ firn^ is obtained. In contrast, the slow-twitch geometry leads to 
a more gradual response, reaching the peak value 3.12 x 10^ firn^ in 43 fis. The "decay 
time" (duration from vesicle release until the decay reaches half of the peak L(/)) for the 
fast-twitch muscle is 213 /xs, significantly shorter than that displayed in the slow-twitch 
geometry, 270 /is. 

In Figure IV.5, A(?) is shown for the slow- and fast-twitch NMJ models. The peak- 
amplitudes are delayed and decay is accelerated in comparison to L(t). As expected, 
the more rapid decay of fast-twitch muscle signal is preserved. There are miniature 
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Figure IV.4: L{t) for NMJ models of fast-twitch and slow-twitch muscles. 

endplate potential (mEPP as opposed to mEPC) data on fast- and slow-twitch muscle in 
mouse [83]. Signal decay according to this data is several times longer, but the paral- 
lels between that geometry and our simulation are incomplete, essentially because the 
respective efforts, in being twenty years apart, are uncoordinated. 

Dystrophic neuromuscular junctions 

Some dystrophic muscle displays improperly developed NMJs in terms of ultrastruc- 
ture [94], but again specific data are limited. For the rectilinear dystrophic NMJ model 
described in Table IV. 1, the rise time, decay time, and amplitude of L(t) are all slightly 
lower than normal (Figure IV.6). Of course, this simply states the obvious: such mus- 
cle does not function as efficiently! Models may not contribute in fundamental ways to 
muscular dystrophy until they can encompass the actual development of NMJ architec- 
ture during generation. 



Reduced AChE density 

As could be said for the distribution of AChR, aberrations in the presence of AChE 
can have functional implications. For the slow-twitch muscle NMJ, lowering AChE 
density to 50 jUm~^ leads to a slightly longer rise time, 46 /xs, with a lower peak level 
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Figure IV.5: A(/) for NMJ models of fast-twitch and slow-twitch muscles. 
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Figure IV.6: Normal versus dystrophic fast-twitch muscle. 



84 




Figure IV.7: Comparing L{t) for different AChE densities. High density: 100 fim ; 



-2. 



—2 



low density: 50 /im ". 



of 2.81 finr' (Figure IV.7). The decay time is about twice as long, at 505 jis. Specific 
experimental information on variations in AChE density among different muscle types 
is limited. As such, this example merely serves to illustrate that an abnormal density 
can have an appreciable effect on the efficiency of synaptic transmission across an NMJ. 

IV.E.3 Realistic neuromuscular junction model 

Finally, we demonstrate the capacity of FEtk to model ACh transmission within a real- 
istic NMJ geometry drawn from electron microscopy. Given this realism, our infrastruc- 
ture is poised to accompany experiments. 



Mesh geometry and initial condition 

We are grateful to have received a realistic 3-dimensional surface-mesh model of an N M J 
from Stiles and coworkers, formed from morphing a 2-dimensional electron micrograph 
[73,74]. The resulting NETGEN volume mesh (Figure IV.8) includes the vesicle, while 
realizing AChE clusters requires the special attention described below. In total, there are 
175 609 vertices and 771 408 simplices in the volume mesh. Edges range from 2.45 x 
10~^ jUm to 1.05 X 10~^ jUm in length, simplex volume ranges from 1.87 x 10~^ jira^ 
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to 2.70 X lO""* jUm^, and the worst (largest) edge ratio in any simplex is 14.4. 

The vesicle is a triangulated sphere of radius 0.024 fim placed roughly in the 
center of the presynaptic membrane. The pore opening is a single triangle, with area 
0.000 54 jUm", at the locus of vesicle attachment. Initially, the ACh concentration is 
zero everywhere but inside the vesicle, where it is 300 mM = 1.8 x 10^ fim~^. This 
translates to approximately 6.9 x 10^ ACh molecules present at the beginning of simu- 
lation. 

AChE clusters and reactivity 

Placing AChE holes into the tetrahedral mesh is a nontrivial process that needs stream- 
lining before "high-throughput" studies can be conducted. Currently, using the APBS 
[95] and FEtk infrastructure, we repeatedly refine, down to a surface area threshold of 
0.002 4 jUm% the simplices that contain the positions of AChE monomers specified by 
Stiles and coworkers. The tetrahedral simplices containing the AChE positions are "re- 
moved" by marking adjacent faces of neighboring simplices as boundary. The resultant 
mesh has 55 589 tetrahedral AChE clusters, distributed as shown in Figure IV.9. The 
average surface area is 0.001 7 /xm^, representative of 59 072 AChE molecules. This 
average cluster surface area must be accounted for when estimating the specific reac- 
tivity k[^^. Since Equation IV.17 scales linearly with surface area of the cluster, the k[^.^ 
estimated in the previous section can be adjusted according to the smaller average sur- 
face area in the realistic model (by a factor of about 1.4). This simplification leads to 
/c^,t = 2.753x10"^ fim-fis-K 

AChR density and postsynaptic detection level 

We use approximately the same AChR density criterion as the above rectilinear mod- 
els. Figure IV.IO shows the postsynaptic detection level L(t) for the more realistic NM J 
model. The rise time is 10 fis, about half as long as that in the previous rectilinear NMJ 
model, but comparison at this stage is somewhat unfounded. 



86 




Figure IV.8: Three views of the realistic mesh: overview, enlarged box containing the 
vesicle on the presynaptic membrane, a cross-section of the secondary cleft (shaded). 
The AChE "holes" are unshaded. The height of the clefts from crest to trough is approx- 
imately 1 jUm. 
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Figure IV.9: Surface area distribution of the AChE holes in the realistic mesh. 
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Figure IV. 10: The postsynaptic detection level L{t) over 150 /is in the realistic NMJ 
model. 
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IV.F Discussion and conclusions 

We have shown that simple rectihnear models are capable of capturing electrophysio- 
logical trends observed in NMJs from different muscle types. Comparing decay time in 
the slow- and fast-twitch muscles reveals a discrepancy in behavior. Also, a model sug- 
gests that dystrophic muscle has a muted postsynaptic response. However, it is likely 
that NMJ malfonnation is less a cause than an effect in the pathology of dystrophy. 
As our primary contribution, we have created a finite element infrastructure capable of 
dealing with a realistic mesh based on electron microscopy data. This motivates future 
partnerships with coordinated experimental investigations on the relationship between 
muscle function and NMJ architecture. Currently, data that correlate geometry and ki- 
netics with muscle function are, when available, scattered and incomplete for our intent 
and purpose. 

Our effort to make physiologically relevant connections is an advancement from 
previous continuum models [71]. The finite element infrastructure we have unveiled 
will evolve on two fronts. First, with improved understanding of reaction kinetics, sim- 
ulations will be more directly related to observable mE P C and mE P P . Second, increasing 
availability of ultrastructural data from electron tomography will offer the capacity for 
extensive comparative studies [96,97]. Only then can parameter space, including NMJ 
geometry, vesicle placement, receptor binding, and reactivity, be comprehensively ex- 
plored to better understand neuromuscular diseases that affect NMJ processes. 
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When a young man succeeds in passing the examinations for even the low- 
est degree, preparations on a most elaborate scale are made at his home for 
honoring him on his return. No one but an eye-witness can imagine the 
scene. A feast is prepared, theatrical performers are often engaged, a pro- 
cession goes out to meet the graduate, who affects all the airs imaginable, 
and his conceit is swollen beyond endurance. His swagger is supercilious to 
the point of silliness. To recognize his old companions is a condescension 
for which they feel extremely grateful. The whole performance tends to 
make these graduates the most obnoxious of all the people one meets. [98] 
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